Avant de commencer…

Contact et échanges

  • Mail :
  • Cours moodle : id 31147, https://moodle.univ-lille.fr/course/view.php?id=31147
    • Mot de passe : gpvnwa
    • Section Données : Les données support pour le cours, les exercices et l’évaluation
    • Section Cours : Présentation et guide mis à jour avant chaque séance. Le support vous aide à suivre le cours, vous pouvez directement copier/coller les codes proposés. Tous les fichiers contenus dans l’archive .zip doivent être extraits dans le dossier du TD sur votre poste pour être lus correctement. Le fichier .html contient l’ensemble de la présentation et s’ouvre via un navigateur web type firefox.
    • Section M1 : Cours de M1 (2021/2022), en particulier pour celles et ceux qui n’étaient pas en M1 ENSP l’an dernier
  • Document partagé (actif jusqu’au 10/04/23) : https://semestriel.framapad.org/p/m2enspr-9wyk, pour la prise de note collaborative, partage de scripts alternatifs, ressources, etc. Vous en faites ce que vous voulez
  • Est-ce que ce mode de fonctionnement vous convient ? Si non, faites une(des) proposition(s) !
  • Avez-vous eu le cours de l’an dernier ? L’avez-vous lu ?

Objectifs du TD

On a vu en M1 ENSP :

  • Apprendre le langage R
  • Importer et mettre en forme des données
  • Présenter des statistiques descriptives et graphiques

Source : allison Horst @allison_horst

Objectifs en M2 ENSP :

  • Assurer les acquis du cours d’initiation de M1
  • Approfondir l’usage du langage R
  • Utiliser des méthodes statistiques plus avancées (/!\ Ce n’est pas un TD de statistiques)
  • Partager des résultats

Source : allison Horst @allison_horst

Organisation

  • 4 demi-journées de 3 heures : 13h-16h15 (dont deux pauses de 15 minutes vers 14h et 15h) les 11, 18 et 25 octobre et 8 novembre en C107 (SH3).
  • 2 demi-journées de 2 heures : 10h-12h (avec une pause de 15 minutes) les 17 et 24 novembre 302 (C1).

Toujours les mêmes bonnes habitudes de travail :

-> Exprimez vos difficultés et blocages : ce qui permet de reprendre ensemble des aspects vus trop rapidement et d’éviter que les difficultés s’accumulent jusqu’à devenir insurmontables

-> Jouez le jeu des exercices : c’est un moyen (le seul) de se familiariser avec l’écriture de R : Les réponses sont sur les supports de cours mais faites l’effort de chercher par vous-même.

Vous (rappels)

  • Que faisiez-vous l’année dernière ? Avez-vous suivi le TD de R - Initiation de M1 ENSP ?
  • Si oui, avez-vous pratiqué depuis le cours de M1 ? Par exemple pour votre mémoire ? Quelles ont été les difficultés ?
  • Si non, avez-vous déjà utilisé ou eu des cours de R ? Dans un autre langage de programmation ?

Programme prévisionnel

  1. Rappels augmentés du M1
  2. Analyses spatiales : cartographie
  3. Données pondérées
  4. Statistiques inférentielles
  5. Production de rapport avec {Rmarkdown}
  6. Selon le temps, création de fonctions avec R

Est-ce que ce programme vous convient et est complémentaire d’autres cours ? Pas de doublons ? D’autres attentes ?

Les différents jeux de données utilisés en cours seront présentés au fur et à mesure.

Evaluation

Il s’agira d’un dossier en groupe de 2 à 4 étudiant·es au choix à partir de données fournies et à rechercher autour d’une même thématique. Le dossier rendu sera :

  • un document (.odt, .pdf ou .html)
  • rédigé de 4 à 6 pages
  • à partir d’un fichier .Rmd
  • comportant au moins une carte
  • présentant des résultats pondérés
  • et au moins un test statistique.

-> Vous m’enverrez le document produit et le script .Rmd au plus tard le 03/01/23 à mon adresse mail (vous aurez un accusé de réception sous 24h).

Barème du dossier

  • Présentation et orthographe /-1 à +1 (ok pour quelques coquilles, à voir comme un bonus)
  • Dossier : /9
    • Présentation des données, de l’enquête /1
    • Problématique, Hypothèses, présentation d’un plan cohérent /2
    • Choix des questions / variables /1
    • Présentation tableaux et graphiques /1
    • Lecture et interprétation tableaux et graphiques /1
    • Développement et analyses en lien avec les questions, données et traitements adéquats pertinence des analyses /1,5
    • Mise en forme {Rmarkdown} /1,5
    • Biblio et références /1
  • Script : /10
    • Importation /1
    • Lisibilité du script, code utile, commentaires explicatifs /1
    • intérêts des traitements /1
    • Données spatiales et cartographie /1
    • Application XXXdes pondérations de la régression /1
    • Finesse et originalité des recodages, mises en forme /1
    • Pertinence des choix graphiques / résultats /1
    • Pas d’erreur de code / logique /1
    • Code efficace /1
    • Bibliothèques, langages hors du cours /1

1. Rappels du TD de M1 - R Initiation

Objectif : Revoir tout le cours de M1 en 3 heures

1.1 Présentation de R et RStudio

R

R est un dérivé gratuit (licence GNU GPL) du langage de statistique S. C’est un projet collaboratif et donc ouvert à tous. Il est entretenu par le R Development Core Team (CRAN). C’est un langage permettant de produire des statistiques descriptives comme des statistiques avancées. Ses fonctionnalités sont étendues ; il offre un large choix d’options graphiques. Il permet de développer ses propres fonctions et d’utiliser des fonctions développées par d’autres chercheurs, ingénieurs ou développeurs.

En bref :

  • Gratuit et opensource
  • Multi-plateformes (Windows, Linux, OS X)
  • Produit simplement des traitements et statistiques avancées :
    • Statistiques descriptives, inférentielles
    • Représentations graphiques
    • Analyses factorielles
    • Modélisations
    • Analyses de réseaux
    • Cartographies
    • Applications, etc.
  • Grâce aux “packages” (extensions, bibliothèques de fonctions), les fonctionnalités sont très étendues et en font un langage complet
  • Importante communauté :
    • Tutoriels et aides sur le web
    • Développement constant de nouvelles fonctionnalités, interfaces et mises à jour régulières
  • Reproductibilité et historique :
    • Partager simplement des rapports commentés (pdf, word, htlm, présentations) et applications en ligne ;
    • Reproduire et adapter des procédures et méthodes réalisées pour d’autres projets (ex : Rapports annuels);
  • Permet de ne pas subir des options par défaut sans les comprendre
  • Compatible avec d’autres langages, exportations graphiques avancées

… Quelques inconvénients :

  • Les fonctionnalités (et bibliothèques) sont très souvent mises à jour et nécessitent d’être maintenues
  • Les fonctionnalités sont créées par les contributeurs et les syntaxes peuvent varier
  • Peine quand les données sont très massives

RStudio (posit)

La console disponible avec R (RGui) est rudimentaire. Il est conseillé d’utiliser un logiciel offrant un espace de travail plus agréable. RStudio est un environnement de développement pour l’écriture et la visualisation de scripts R. Le logiciel est produit par une société (qui va prochainement s’appeler posit) à l’origine de nombreuses ressources et projets autour de R ({Rmarkdown} pour les rapports, {shiny} pour les applications, {tidyverse}, RstudioCloud, etc.). Il est également disponible gratuitement. Il offre un éditeur de texte (coloration, auto complétion, indentation) et un affichage simultané des scripts, sorties, fichiers, graphiques, aides. Il est quasiment systématiquement utilisé pour écrire en R.

Téléchargements

- (1) Télécharger le langage R :

- (2) Télécharger et installer l’IDE RStudio :

Environnement de RStudio

Attention :

  • Respect de la casse
  • Ni accent, ni espace, ni caractère spécial en-dehors de “_”

Les packages

  • Packages : Extensions ajoutées à votre session de travail, une fois installées et chargées. Les packages sont des bibliothèques de fonctions supplémentaires développées par d’autres personnes.

  • Les packages “officiels” sont validés par la communauté de R (CRAN). Ils sont tous accompagnés d’une documentation (document pdf issu du site cran.r-project.org).

Pour utiliser un package validé par le CRAN :

  • 1. L’installer (une seule fois)
install.packages("tidyverse", dep=T)

Cette étape d’installation des packages est délicate : R évoluant rapidement, la compatibilité entre les packages et entre les packages et la version utilisée de R n’est toujours maintenue.

  • 2. Le charger (à refaire pour chaque session).
  • En ligne de commande :
library(tidyverse)
library(questionr)

/!\ Il est préférable d’écrire ces lignes afin que les packages soient automatiquement définis et que les fonctions associées fonctionnent.

On recommande aussi de mettre ces lignes en début de script.

library() est équivalent à require()

Enregistrer votre script : File > Save as…

Organiser son travail

Quand on entame un projet, les lignes de commande, scripts, tableaux créés et dossiers peuvent s’accumuler. Il est important d’organiser son travail dès le début.

  • Gérer ses projets. Un projet est une session qu’on associe à un dossier :

    • Vous retrouvez votre session : script, objets, recodages et transformations déjà effectués

    • Renvoie directement à cet emplacement

    • Evite de mélanger les données et scripts portant sur différents sujets

    • Permet de partager l’ensemble du dossier à d’autres personnes

En haut à droite, une icône en forme de cube bleu permet de créer un “nouveau projet”. C’est un assistant de création d’un dossier dédié à un projet spécifique dans votre explorateur (windows, mac ou linux). Il sera associé à une session de Rstudio, qui sera indépendante des documents et objets utilisés pour un autre projet.

ou bien File > New Project > Existing directory –> Browse > Sélectionnez votre dossier créé pour le TD > Create project

Un fichier Nom_du_dossier.Rproj apparaît dans le dossier, c’est le raccourci qui mènera à la session R. Un fichier .RData sera également créé et contiendra les tableaux présents dans l’environnement.

  • Commenter ses scripts

    • Les commentaires sont précédés de # et ne sont pas exécutés.

    • Pour passer un ensemble sélectionné en commentaire (et inversement) : ctrl+shift+c.

    • Aérer le texte et mettre des titres facilitent aussi grandement la lecture de scripts : sur une ligne, # Titre ----

  • Dernier rappel

Il existe toujours plusieurs moyens d’aboutir à un résultat ; L’essentiel est de comprendre ce que vous recopiez / adaptez / écrivez et de favoriser les solutions les plus simples (pour faciliter la relecture, minimiser la mémoire sollicitée et les risques d’erreurs)

  • Veillez à :

    • Organiser vos différents scripts à l’aide des sections et commentaires : un·e futur·e vous ou un·e collègue doit pouvoir se replonger facilement dans le code

    • Ne laisser que le code utile ne comportant pas d’erreur ou d’instructions longues à exécuter et/ou inutile à exécuter (comme un install.packages). Alternez entre les zones console (retours et tests) et les script (code enregistré). Dans l’idéal, vous devriez pouvoir exécuter l’ensemble du script sans erreur (le bouton “source” exécute l’ensemble du script et s’arrête à la première erreur).

-> Si vous ne l’avez pas fait, créez un dossier consacré au TD de R avancé (M2), appelé par exemple “TD_R_M2”.

-> Téléchargez le support de cours pour le mettre dans ce dossier.

-> Créez un projet .Rproj associé au dossier du TD et ouvrez un script .R

Fichier > Nouveau projet > Existing directory > Browse > Choisir le dossier créé pour le TD

1.2. Exercices avec les données Coco1

1. Consultation de la documentation

Avant toutes choses, commencez à lire la documentation associée aux premières données Coco1 qu’on utilisera dans les parties 1 à 3 du TD :

  • Documentation_coco1.pdf : pp. 4-8 ; La suite est le dictionnaire des codes

    • Quels apports et originalités ?
    • Quelles zones d’ombre et limites ?

  • Questionnaire_coco1.pdf : L’ensemble ; En notant les codes des questions qui vous intéressent

  • !! Attention : Distinguez bien la Documentation qui comporte le dictionnaire des codes du Questionaire qui ne contient que les questions de l’enquête. Les données reçues contiennent en plus des réponses au questionnaire une série de variables socio-démographiques issues de l’enquête annuelle réalisée par l’équipe ELIPSS. Vous avez toujours le fichier Documentation_coco1.pdf ouvert pendant que vous travaillez sur les données Vous trouverez le questionnaire de l’enquête annuelle dans l’archive cdsp_ea2019, fichier Questionnaire_EA2019.pdf.

2. Importation des données coco de la vague 1

Rappels - Importation des données

Suivant les types des données à importer (voir leur extension), différentes fonctions d’importation seront utilisées :

  • Pour les fichiers .rds : readRDS("chemin/nomFichier.rds") ;
  • Pour les fichiers texte .csv : read.csv2("chemin/nomFichier.csv", stringsAsFactors=F) ;
  • Pour les fichiers .txt : read.table("chemin/nomFichier.txt", sep="\t")

Un package utile facilitant les importations et exportations (entre la session R et le dossier associé au projet) : rio : voir la vignette de présentation

… /!\ : Beaucoup de packages facilitent les actions : risque de s’y perdre et risque d’être facilement bloqué·e en cas de bug si vous ne connaissez pas un minimum le R-base

unzip("data/cdsp_coco1_2020/cdsp_coco1_2020.zip", exdir="data/")
coco <- read.csv("data/cdsp_coco1_2020/donnees_coco1/csv_coco1/fr.cdsp.ddi.elipss.202004coco1.csv",
                 stringsAsFactors = F)

- Repérez les noms des variables contenant le genre, l’âge et le niveau de diplôme des enquêté·es.

Le genre se trouve dans 2 variables (eayy_A1 et coco1_cal_SEXE), l’âge dans plusieurs (eayy_A2A_rec pour l’âge quinquennal et coco1_cal_AGE1 pour l’âge décennal) et le niveau de diplôme est dans la variable eayy_B18_rec.

La différence est que les variables préfixées eayy sont les variables issues de l’enquête annuelle (à utiliser) et les variables préfixées cal sont les variables de calages utilisées pour produire les coefficients de redressement et qu’on utilisera pas. On utilisera plutôt les variables eayy.

  • On remarquera que les codes ne sont pas des plus explicites dans cette table de données.
  • Il faudra lire en complément du questionnaire le détail des variables socio-démographiques issues de l’enquête annuelle et qui se trouve dans la Documentation_coco1.pdf pp.83-130.

Rappels - Calculs sur les tables

  • Effectifs avec table :

    • Tris à plat : table(df$Variable)
    • Tris croisés : table(df$1èreVariable, df$2èmeVariable)

  • Effectifs et poucentages sur une variable avec freq : freq(df$variable)

  • Pourcentages avec prop.table(), rprop() et cprop() s’appliquant sur une table d’effectifs : rprop(table(df$1èreVariable, df$2èmeVariable))

  • Indices de répartition : min, max, sd, mean, median et summary.

  • Exporter un résultat : write.csv2(resultat, "resultats/tc_prepa_genre.csv") ; s’assurer qu’il existe un dossier “resultat” puis mise en forme sous Excel ou Calc

3. Quelle part de femmes et d’hommes ont répondu au questionnaire ? Présentez un tableau clair et lisible

Rappels - Recodages

Recodages par assignation :

table$sexe[table$sexe=="1"] <- "Homme"

Recoder avec {tidyverse} :

# 2. Recodage :
table$sexe <- fct_recode(
  factor(table$sexe),
  "Homme"="1")
# ou fct_collapse()

Recoder en NA et inversement :

table$recodage[table$recodage=="Non renseigné"]<-NA
table$recodage[is.na(table$recodage)]<-"Non renseigné"
# Ou
table$recodage <- fct_recode(factor(table$variable_a_recoder),
                             "Valeur 1"="1","Valeur 2"="2",
                     NULL="Non renseigné")
table$recodage <- fct_explicit_na(table$recodage,
                                 na_level = "Non renseigné")

# 1. Chargement des packages utilisés :
library(questionr)
library(tidyverse)
# 2. Recodage :
coco$sexe <- fct_recode(
  factor(coco$eayy_A1),
  "Homme"="1", "Femme"="2")
# 3. Tableau de fréquences :
freq(coco$sexe,
     total= T, valid=F)

Recodages conditionnels : créer de nouvelles variables à partir de celles existantes

2 fonctions principales permettent de recoder des variables sur conditions :

ifelse(condition, si condition remplie, si condition non remplie)
case_when( conditions 1 ~ effet de la condition 1 ,
condition 2 ~ effet de la condition 2,
etc,
TRUE ~ effet si aucune condition n'est remplie)

4. Faire le tableau et le graphique associé du sentiment de tristesse selon le genre des enquêté·es

# 1. Recodages :
coco$coco1_q19 <- fct_recode(factor(coco$coco1_q19),
                             "En permanence"="1",
                             "Souvent"="2",
                             "Quelques fois"="3",
                             "Rarement"="4",
                             "Jamais"="5",
                             NULL = "9999")
# 2. Tableau :
tc_genre_tristesse <- rprop(table(sexe=coco$sexe,
            tristesse=coco$coco1_q19))
tc_genre_tristesse
##           tristesse
## sexe       En permanence Souvent Quelques fois Rarement Jamais Total
##   Homme      0.8           3.8    18.0          36.5     40.9  100.0
##   Femme      1.1           9.5    32.8          29.6     27.0  100.0
##   Ensemble   0.9           6.8    25.9          32.9     33.5  100.0

Rappel - Graphiques avec {ggplot2}

Source : allison Horst @allison_horst

GGplot2 est un package dédié aux graphiques. Sa syntaxe est cohérente avec celle de Tidyverse. L’intérêt de Ggplot2 est de produire rapidement et efficacement des graphiques complexes par ajout de couches successives. De plus, il est possible d’enchaîner la mise en forme des tableaux avec la réalisation des graphiques.

La syntaxe de base de ggplot est :

ggplot (table_de_données) + aes (x = Variable_abscisses, y = Variable_ordonnées ) + geom_X (elements de mise en forme) + labs (titres)

où ggplot() définit la source des données, geom_X() le type de représentation, aes() les données à représenter et labs() les différents titres et textes.

A laquelle on peut ajouter des précisions du type :

  • theme() : Mise en forme des zones de texte et de dessin

  • scale_y_continuous() : Options sur les échelles

  • scale_fill_manual() : Options sur les couleurs

# 4. Graphique :
tc_genre_tristesse %>% 
  data.frame() %>% 
  filter(tristesse !="Total") %>% 
  ggplot()+
  aes(fill=sexe, x=tristesse, y= Freq) %>% 
  geom_bar(col="black",stat="identity", position="dodge")+
  labs(title="Tristesse selon le genre pendant le 1er confinement",
       x="S'est senti·e triste",y="%",fill="",
       caption = "Source : Coco1 - ELIPSS, CDSP, avril 2020")+
  viridis::scale_fill_viridis(discrete=T)+
  theme_bw()

5. Quelle est la part moyenne des personnes de l’entourage des enquêté·es qui ont eu le covid et sont aujourd’hui guéries ? Et la part médiane ?

On utilise les questions Q12 et Q13. Il y a peu de variables réellement numériques.

class(coco$coco1_q12)
## [1] "integer"
coco$coco1_q12[coco$coco1_q12 %in% c(6666,9999)] <- NA
coco$coco1_q13[coco$coco1_q13 %in% c(6666,9999)] <- NA
coco$gueries <- round(coco$coco1_q13/coco$coco1_q12*100,1)
summary(coco$gueries)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##    0.00   33.30  100.00   65.52  100.00  100.00     648

En fonction de l’âge ? En fonction de la PCS ?

Rappels - Mise en forme de tables avec {tidyverse} et %>%

Sélections, filtres et tris

Les principales fonctions permettant de réarranger les tables via tidyverse sont :

  • select(table, colonnes à conserver) : Sélection de colonnes ; Voir starts_with(), contains(), ends_with()
  • filter(table, condition du filtre) : Sélection de lignes
  • arrange(table, colonnes sur laquelle faire le tri, desc(2ème colonne sur laquelle faire un tri décroissant)) : Tri de la table
  • slice(table, position des lignes à conserver) : Sélection de lignes
  • mutate(table, nouvelleVariable = ) : Création de variables
  • rename("nouveau nom" = "ancien nom") : Renommer une variable
  • L’opérateur %>% (pipe, ctrl+shift+m) permet d’enchaîner des fonctions : f(a) s’écrit a %>% f().

Rappels - Regroupements de lignes

Il peut être intéressant de regrouper plusieurs observations selon un ou plusieurs critères afin de produire des calculs. La fonction group_by permet de produire ces regroupements et est associée à une fonction de calcul dans summarise.

Exemple : On souhaite compter le nombre d’individus par classe d’âge ainsi que le revenu moyen par classe d’âge d’une table donnée :

table %>% 
  group_by(classe_dage) %>% 
  summarise(Nb_par_classe_dage = n(),
            Revenu_moyen = mean(revenu, na.rm=T))

En fonction de l’âge ?

# 1. Recodage :
coco$age_quinq <- fct_recode(factor(coco$eayy_A2A_rec),
                             "<24 ans"="4", "25-29 ans"="5",
                             "30-34 ans"="6", "35-39 ans"="7",
                             "40-44 ans"="8", "45-49 ans"="9",
                             "50-54 ans"="10", "55-59 ans"="11",
                             "60-64 ans"="12","65-69 ans"="13",
                             "70 ans et+"="14")
# 2. Avec group_by :
coco %>% 
  group_by(age_quinq) %>% 
  summarise(mediane = median(gueries, na.rm=T),
            moyenne = mean(gueries, na.rm=T))

En fonction de la PCS ?

coco$pcs <- fct_collapse(factor(coco$eayy_PCS6CJT),
                            "Agriculteurs"="1",
                       "Artisans, comm., chef d'E"="2",
                       "Cadres et pr.intell.sup."="3",
                       "Prof. intermédiaires"="5",
                       "Employés"="4","Ouvriers"="6",
                       NULL = c("666","999"))

coco %>% 
  group_by(pcs) %>% 
  summarise(mediane = median(gueries, na.rm=T),
            moyenne = mean(gueries, na.rm=T))

6. Retirez les préfixes “coco1_” des noms de variables de l’enquête

names(table) permet d’agir sur les noms des variables

gsub(pattern = "on", replacement = "a", x = "bonbon") modifie tous les “on” en “a” de “bonbon”.

names(coco) <- gsub("coco1_","",names(coco))

7. On veut maintenant comprendre l’effet du 1er confinement (et du début de la crise sanitaire) sur la situation d’emploi des enquêté·es

On passera par la création de variables intermédiaires :

  • sitpro.avt : Situation d’emploi avant le confinement ; A partir de la variable q32

  • sitpro.conf : Situation d’emploi pendant le confinement ; A partir de la variable q37

  • tt.conf : Et on intègre le télétravail avec q38

sitpro.avt

# freq(coco$q32)
coco$sitpro.avt <- fct_collapse(factor(coco$q32),
                            "En emploi, contrat" = c("1","2"),
                            "Sans emploi" = c("3","4"),
                            "Retraite" = "6",
                            "Autre inactif·ve" = c("5","7","8"), NULL="9999")
freq(coco$sitpro.avt, sort="dec", valid=F)

sitpro.conf et tt.conf

# freq(coco$q37)
coco$sitpro.conf <- fct_collapse(factor(coco$q37),
                                        "En congé"=c("1","2","3"),
                                        "Au chômage"=c("4","5"),
                                        "En emploi"=c("6"),
                                 NULL=c("7","6666","9999"))
freq(coco$sitpro.conf)

sitpro.conf et tt.conf

# freq(coco$q38)
coco$tt.conf <- fct_recode(factor(coco$q38),
                                 "Lieu habituel"="1",
                                 "Lieu habituel/TT"="2",
                                 "Télétravail"="3", NULL = "6666")
freq(coco$tt.conf)

sitpro.tt.conf : Situation d’emploi et TT dans une colonne

# table(coco$sitpro.conf, coco$tt.conf, useNA = "ifany")
coco$sitpro.tt.conf <- ifelse(is.na(coco$tt.conf) ==F,
                           paste0("En emploi - ", coco$tt.conf),
                           as.character(coco$sitpro.conf))
freq(coco$sitpro.tt.conf)

De fait, les questions sur leschangements de situation d’emploi ne concernent que les personnes précédemment en emploi. Pour analyser ces questions, on passe par la création d’une table conservant la sous-population des actif·ves avant la pandémie.

–

coco.activite <- coco %>% filter(sitpro.avt == "En emploi, contrat")
# freq(coco.activite$sitpro.tt.conf, sort="dec",valid=F)

On peut maintenant envisager les modifications des conditions de travail selon l’âge, la PCS, etc.

rprop(table(coco.activite$pcs, coco.activite$sitpro.tt.conf,
            useNA="ifany"))
##                            
##                             Au chômage En congé En emploi
##   Agriculteurs               25.0       12.5      0.0    
##   Artisans, comm., chef d'E  21.1       15.8      0.0    
##   Cadres et pr.intell.sup.   17.7        6.3      2.5    
##   Employés                   18.0       15.3      7.2    
##   Prof. intermédiaires       24.8       20.0      3.8    
##   Ouvriers                   29.0       14.5      4.3    
##   <NA>                       25.6       10.6      6.0    
##   Ensemble                   23.2       13.2      4.9    
##                            
##                             En emploi - Lieu habituel
##   Agriculteurs               50.0                    
##   Artisans, comm., chef d'E  31.6                    
##   Cadres et pr.intell.sup.   20.3                    
##   Employés                   19.8                    
##   Prof. intermédiaires       29.5                    
##   Ouvriers                   30.4                    
##   <NA>                       21.6                    
##   Ensemble                   24.2                    
##                            
##                             En emploi - Lieu habituel/TT
##   Agriculteurs                0.0                       
##   Artisans, comm., chef d'E  10.5                       
##   Cadres et pr.intell.sup.    6.3                       
##   Employés                    5.4                       
##   Prof. intermédiaires        1.9                       
##   Ouvriers                    2.9                       
##   <NA>                        5.5                       
##   Ensemble                    4.7                       
##                            
##                             En emploi - Télétravail <NA>  Total
##   Agriculteurs                0.0                    12.5 100.0
##   Artisans, comm., chef d'E  21.1                     0.0 100.0
##   Cadres et pr.intell.sup.   45.6                     1.3 100.0
##   Employés                   32.4                     1.8 100.0
##   Prof. intermédiaires       17.1                     2.9 100.0
##   Ouvriers                   17.4                     1.4 100.0
##   <NA>                       29.1                     1.5 100.0
##   Ensemble                   27.8                     1.9 100.0

!! Les effectifs commencent à devenir faibles (en lien avec la structure de la population du panel ELIPSS).

freq(coco.activite$pcs)

Des catégories devraient être regroupées…

8. Risques professionnels des personnes en activité selon le genre

Pour en terminer avec les exercices de rappel du M1, on veut représenter graphiquement les risques professionnels encourus pendant la période de confinement pour les personnes déjà en activité en fonction du genre.

Les variables vont les q41_1 à q41_6 et viennent d’une question à choix multiples :

Si Q38 = {1, 2}

Q41 - Au cours des deux dernières semaines, vos conditions de travail vous ont-elles exposé à : (Plusieurs réponses possibles)

  1. Des risques pour ma santé
  2. Des conflits avec mon employeur
  3. Des conflits avec mes collègues
  4. Un risque de chômage
  5. Une perte de salaire
  6. Aucune de ces situations [Exclusif] [EMPTY+1]

.. plutôt dommage qu’on ne pose pas la question aux personnes en télétravail.

On commence par les recodages :

# Recodage du risque pour la santé :
coco.activite$risque.sante <- fct_recode(factor(coco.activite$q41_1),
                              "Oui"="1","Non"="2",NULL = "6666")
freq(coco.activite$risque.sante)
# Recodage du risque de conflit avec l'employeur :
coco.activite$risque.conflit.empl <- fct_recode(factor(coco.activite$q41_2),
                              "Oui"="1","Non"="2",NULL = "6666")
freq(coco.activite$risque.conflit.empl)

Poursuivez les recodages…

# Recodages suivants :
coco.activite$risque.conflit.coll <- fct_recode(factor(coco.activite$q41_3),
                              "Oui"="1","Non"="2",NULL = "6666")
coco.activite$risque.cho <- fct_recode(factor(coco.activite$q41_4),
                              "Oui"="1","Non"="2",NULL = "6666")
coco.activite$risque.baisse.sal <- fct_recode(factor(coco.activite$q41_5),
                              "Oui"="1","Non"="2",NULL = "6666")
coco.activite$risque.aucun <- fct_recode(factor(coco.activite$q41_6),
                              "Oui"="1","Non"="2",NULL = "6666")

Représentez une de ces variables en graphique :

coco.activite %>% 
  filter(is.na(risque.sante)==F) %>% 
ggplot()+
  aes(fill=risque.sante, x=sexe)+
  geom_bar(position="fill")

Ensuite ?

On essayera de présenter tous ces croisements “collés” dans un unique tableau.

rprop(table(sexe=coco.activite$sexe, risque=coco.activite$risque.sante))
##           risque
## sexe       Oui   Non   Total
##   Homme     47.8  52.2 100.0
##   Femme     63.7  36.3 100.0
##   Ensemble  57.3  42.7 100.0

On utilise rbind qui “colle” 2 tables en lignes. Avec 2 risques :

risques <-  rbind(
   # 1er tableau croisé :
rprop(table(sexe=coco.activite$sexe, risque=coco.activite$risque.sante)) %>% 
  data.frame() %>% 
  mutate(type="Santé"),
# 2ème tableau croisé :
 rprop(table(sexe=coco.activite$sexe, risque=coco.activite$risque.conflit.empl)) %>% 
  data.frame() %>% 
  mutate(type="Conflit employeur"))

Puis mise en forme et graphique :

risques %>% 
  filter(risque=="Oui" & sexe !="Ensemble") %>% 
  ggplot()+
  aes(x=type, y=Freq)+
  geom_bar(stat="identity")+
  facet_wrap(~sexe)

Poursuivre avec l’ensemble des risques

risques <- rbind(risques,
                rprop(table(sexe=coco.activite$sexe, risque=coco.activite$risque.conflit.coll)) %>% 
                  data.frame() %>% 
  mutate(type="Conflit collègues"))
risques <- rbind(risques,
                rprop(table(sexe=coco.activite$sexe, risque=coco.activite$risque.cho)) %>% 
  data.frame() %>% 
  mutate(type="Chômage"))
risques <- rbind(risques,
                rprop(table(sexe=coco.activite$sexe, risque=coco.activite$risque.baisse.sal)) %>% 
  data.frame() %>% 
  mutate(type="Baisse de salaire"))
risques <- rbind(risques,
                rprop(table(sexe=coco.activite$sexe, risque=coco.activite$risque.aucun)) %>% 
  data.frame() %>% 
  mutate(type="Aucun"))

Et le graphique :

risques %>% 
  filter(risque=="Oui" & sexe !="Ensemble") %>% 
  ggplot()+
  aes(x=type, fill=sexe, y=Freq)+
  geom_bar(stat="identity",position="dodge")+
  theme(axis.text.x = element_text(angle=45, hjust=1))+
  labs(title="Exposition aux risques professionnels", fill="Genre",
       x="",y="%", 
       caption = "Champ : Personnes en activité\nSource : Coco1 - ELIPSS, CDSP, avril 2020")

Et la table :

risques %>% 
  filter(risque =="Oui") %>% 
  mutate(Freq = round(Freq,1)) %>% 
  pivot_wider(-risque,names_from = "sexe", values_from ="Freq")

Rappels sur les jointures de tables

On l’a également vu l’an dernier mais un rapide rappel de la façon de joindre différentes tables vous aidera pour le dossier d’évaluation.

L’enquête “Faire face au covid19” s’est faite en plusieurs vagues. Tou·tes les enquêté·es n’ont pas répondu à chaque vague, il faudra donner certaines précisions sur le type de jointure attendu.

coco2 <- read.csv("data/cdsp_coco2_2020/donnees_coco2/csv_coco2/fr.cdsp.ddi.elipss.202004.csv",
                   stringsAsFactors = F)
jointure.left <- left_join(coco, coco2,
                    by = "UID_ANONYM_COVID")
jointure.right <- right_join(coco, coco2,
                    by = "UID_ANONYM_COVID")
jointure.full <- full_join(coco, coco2,
                    by = "UID_ANONYM_COVID")
jointure.inner <- inner_join(coco, coco2,
                    by = "UID_ANONYM_COVID")
nrow(jointure.left)
## [1] 1076
nrow(jointure.right)
## [1] 998
nrow(jointure.full)
## [1] 1166
nrow(jointure.inner)
## [1] 908

Poursuivre les exercices de rappels du M1

D’ici la prochaine séance, si vous souhaitez poursuivre les exercices, vous pouvez reprendre le déroulé du corrigé de l’examen de l’an dernier (que je vous avais déjà envoyé l’an dernier en M1) -> Section “M1 R Initiation (2020/2021)” > Examen. On pourra y revenir rapidement en début de prochaine séance si vous avez des questions.

1.3. Retour sur la correction de l’examen de M1 de 2021/2022

1. Créez un dossier et un projet R pour l’évaluation, créez un script R, chargez le package « tidyverse » et importez les 4 tables en lignes de commande.

# Chargement du package :
library(tidyverse)
# Importation des données :
etab <- read.csv2("data/Examen_M1_21-22/fr-esr-principaux-etablissements-enseignement-superieur.csv",
                  stringsAsFactors = F, fileEncoding = "UTF-8")
etu <- read.csv2(
  "data/Examen_M1_21-22/fr-esr-statistiques-sur-les-effectifs-d-etudiants-inscrits-par-etablissement.csv",
                  stringsAsFactors = F, fileEncoding = "UTF-8")
perm <- read.csv2("data/Examen_M1_21-22/fr-esr-enseignants-titulaires-esr-public.csv",
                  stringsAsFactors = F, fileEncoding = "UTF-8")
nperm <- read.csv2("data/Examen_M1_21-22/fr-esr-enseignants-nonpermanents-esr-public.csv",
                  stringsAsFactors = F, fileEncoding = "UTF-8")

La principale difficulté a été la non-prise en compte de l’encodage. !! Vérifiez toujours que les données ont correctement été importées. Les aspects à prendre en compte :

  • Le format du fichier :

-> Conditionne la fonction à utiliser pour lire le contenu du fichier et le stocker dans un objet

  • Le caractère qui séparera les cellules :

-> Préciser si besoin le caractère séparateur et Vérifier avec names(table) ou str(table)

  • L’encodage des caractères :

-> En version numérique, les caractères sont stockés sous forme de codes en octets. Notamment en raison de différents alphabets, il existe plusieurs systèmes de reconnaissance de ces codes : C’est l’encodage. En plus de ça, les systèmes d’exploitation ont des systèmes d’encodage par défaut différents : Linux / Mac <-> UTF-8 ; Windows (ordinateur configuré en langue Française) <-> WINDOWS-1252 (ou ISO-8859-1). UTF-8 est plus universel et est plus standard. Beaucoup de données sont fournies avec ce système d’encodage.

-> Préciser le système d’encodage et vérifier que les caractères spéciaux (accentués, par exemple) sont correctement lisibles avec unique(table$VariableNominale) par exemple.

Table Établissements de l’enseignement supérieur en France :

2. En utilisant la table « Principaux établissements d’enseignement supérieur », donnez le nombre moyen d’étudiant·es dans l’ensemble des établissements en 2017/2018.

Quelle était la difficulté ?

class(etab$Effectifs.d.étudiants.inscrits.2017.18)
## [1] "character"
head(etab$Effectifs.d.étudiants.inscrits.2017.18)
## [1] "2 333" "1 101" ""      "986"   "6 103" "2 154"

La colonne contenant le nombre d’étudiant·es dans les établissements n’était pas de type numérique. De plus, elle comportait des espaces insécables séparant les milliers. Il arrive que les données obtenues ne soient pas formatées comme on l’aurait souhaité. Comment la reformater correctement ?

etab$etu_1718 <- gsub(" ","",etab$Effectifs.d.étudiants.inscrits.2017.18)
etab$etu_1718 <- as.integer(etab$etu_1718)
mean(etab$etu_1718, na.rm=T)
## [1] 8593.155

3. Quel est le nombre moyen d’étudiant·es en 2017/2018 selon le secteur (privé ou public) des établissements ?

En version R-base :

tapply(etab$etu_1718, etab$secteur_d_etablissement, mean, na.rm=T)
##               Privé    Public 
##       NaN  2686.532 12047.972

En version tidyverse :

etab %>% 
  group_by(Secteur=secteur_d_etablissement) %>% 
  summarise(`Nb moyen d'étudiant·es en 17/18` = mean(etu_1718, na.rm=T))

On aurait aussi pu ajouter le nombre d’établissements pour chacun des secteurs et arrondir la moyenne :

etab %>% 
  group_by(Secteur=secteur_d_etablissement) %>% 
  summarise(`Nb moyen d'étudiant·es en 17/18` = round(mean(etu_1718, na.rm=T),0),
            "Nombre d'établissements"=n())

4. Créez une variable, toujours dans la table établissements, distinguant les établissements par tranches de date de création (Avant 1900, 1900-1944, 1944-1969, 1970-1999, 2000 ou après).

etab$date_creation_cl <- case_when(is.na(etab$date_creation) ~ "NR",
                                   etab$date_creation < 1900 ~ "Avant 1900",
                                   etab$date_creation < 1945 ~ "1900-1944",
                                   etab$date_creation < 1970 ~ "1945-1969",
                                   etab$date_creation < 2000 ~ "1970-1999",
                                   T ~ "2000 ou après")
# Mettre « Avant 1900 » en 1er :
etab$date_creation_cl <- fct_relevel(etab$date_creation_cl, "Avant 1900")
library(questionr)
freq(etab$date_creation_cl)

5. Quelle part d’établissements dispose d’un compte Youtube ?

freq(etab$compte_youtube!="")

6. Où trouve-t-on le plus d’établissements privés ?

addmargins(table(etab$Ancienne.région, etab$secteur_d_etablissement))
##                             
##                                  Privé Public Sum
##                                0     0      6   6
##   Alsace                       0     0      7   7
##   Aquitaine                    0     2      6   8
##   Auvergne                     0     0      2   2
##   Basse-Normandie              0     1      3   4
##   Bourgogne                    0     1      2   3
##   Bretagne                     0     5     13  18
##   Centre-Val de Loire          0     0      3   3
##   Champagne-Ardenne            0     3      3   6
##   Corse                        0     0      1   1
##   Franche-Comté                0     1      3   4
##   Guadeloupe                   0     0      1   1
##   Guyane                       0     0      1   1
##   Haute-Normandie              0     3      6   9
##   Île-de-France                1    35     53  89
##   La Réunion                   0     0      1   1
##   Languedoc-Roussillon         0     1      4   5
##   Limousin                     0     1      1   2
##   Lorraine                     0     2      1   3
##   Mayotte                      0     0      1   1
##   Midi-Pyrénées                0     3     16  19
##   Nord-Pas-de-Calais           0     7      5  12
##   Pays de la Loire             0    11      4  15
##   Picardie                     0     4      3   7
##   Poitou-Charentes             0     2      3   5
##   Provence-Alpes-Côte d'Azur   0     2      7   9
##   Rhône-Alpes                  0     5     14  19
##   Sum                          1    89    170 260
rprop(table(etab$Ancienne.région, etab$secteur_d_etablissement))
##                             
##                                    Privé Public Total
##                                0.0   0.0 100.0  100.0
##   Alsace                       0.0   0.0 100.0  100.0
##   Aquitaine                    0.0  25.0  75.0  100.0
##   Auvergne                     0.0   0.0 100.0  100.0
##   Basse-Normandie              0.0  25.0  75.0  100.0
##   Bourgogne                    0.0  33.3  66.7  100.0
##   Bretagne                     0.0  27.8  72.2  100.0
##   Centre-Val de Loire          0.0   0.0 100.0  100.0
##   Champagne-Ardenne            0.0  50.0  50.0  100.0
##   Corse                        0.0   0.0 100.0  100.0
##   Franche-Comté                0.0  25.0  75.0  100.0
##   Guadeloupe                   0.0   0.0 100.0  100.0
##   Guyane                       0.0   0.0 100.0  100.0
##   Haute-Normandie              0.0  33.3  66.7  100.0
##   Île-de-France                1.1  39.3  59.6  100.0
##   La Réunion                   0.0   0.0 100.0  100.0
##   Languedoc-Roussillon         0.0  20.0  80.0  100.0
##   Limousin                     0.0  50.0  50.0  100.0
##   Lorraine                     0.0  66.7  33.3  100.0
##   Mayotte                      0.0   0.0 100.0  100.0
##   Midi-Pyrénées                0.0  15.8  84.2  100.0
##   Nord-Pas-de-Calais           0.0  58.3  41.7  100.0
##   Pays de la Loire             0.0  73.3  26.7  100.0
##   Picardie                     0.0  57.1  42.9  100.0
##   Poitou-Charentes             0.0  40.0  60.0  100.0
##   Provence-Alpes-Côte d'Azur   0.0  22.2  77.8  100.0
##   Rhône-Alpes                  0.0  26.3  73.7  100.0
##   Ensemble                     0.4  34.2  65.4  100.0
# En nombres en Île-de-France – où se trouvent 1/3 des établissements
# de France (freq(etab$Ancienne.région=="Île-de-France")) -, en proportion
# dans les Pays de la Loire.

7. Sortez dans la console les noms, secteurs (privé ou public) et le nombre d’étudiant·es des 5 établissements comptant le plus d’étudiant·es (en 17/18) dans le Nord-Pas-de-Calais. Que remarquez-vous ?

etab %>% 
  filter(Ancienne.région =="Nord-Pas-de-Calais") %>% 
  arrange(desc(etu_1718)) %>% 
  select(Nom=Libellé, Secteur=secteur_d_etablissement, 
         `Nb d'étudiant·es`=etu_1718, `Date de création`=date_creation) %>% 
  slice(1:5)

Table des Effectifs étudiants

8. Ne conservez dans cette table que les données relatives à l’année 2020/2021.

etu <- etu %>% filter(Année.universitaire=="2020-21")

9. Ajoutez pour chaque établissement une colonne présentant la part de femmes parmi les étudiant·es et une colonne comprenant la part d’étudiant·es en sciences de l’ingénierie.

# names(etu)
etu <- etu %>% 
  rename(etu_total = 
Nombre.d.étudiants.inscrits..inscriptions.principales..y.compris.doubles.inscriptions.CPGE) %>% 
  mutate(ratioFemmes = round(Sexe...Femmes/etu_total*100,2),
         ratioInge = 
round(Grande.discipline...Sciences.et.sciences.de.l.ingénieur/etu_total*100,2))

10. Faites un graphique présentant au mieux ces deux données (la part de femmes et la part d’étudiant·es en Sciences et sciences de l’ingénierie).

etu %>% 
  filter(ratioInge >0) %>% 
  ggplot() +
  aes(x = ratioFemmes, y = ratioInge, size = etu_total,
      col = Type.d.établissement) +
  geom_point(alpha=.6) +
  theme_bw()+
  labs(title = "Femmes en sciences de l'ingénierie - 20/21",
       x = "Part de femmes (%)", y = "Part en sciences de l'ingénierie (%)",
       caption = "Source : Ministère de l'ESRI, 2021",
       size = "Nombre d'étudiant·es",
       col = "Type d'étab.")

11. Dans un tableau, présentez le nombre d’enseignant·es permanent·es, non-permanent·es et la part de non-permanent·es dans les établissements des Hauts-de-France en 2019/2020. Attention, c’est un tableau de données agrégées, les effectifs sont dans la colonne « effectif ».

full_join(perm %>% 
  filter(Rentrée == "2019") %>% 
  filter(Région=="Hauts-de-France") %>% 
  group_by( Établissement) %>% 
  summarise(permanents=sum(effectif)),
nperm %>% 
  filter(Rentrée == "2019") %>% 
  filter(Région=="Hauts-de-France") %>% 
  group_by( Établissement) %>% 
  summarise(`non-permanents`=sum(Effectif)),
by="Établissement") %>% 
  mutate(`Part de non-permanent·es` = 
      round(`non-permanents`/(`non-permanents`+permanents)*100,2)) %>% 
  arrange(desc(`Part de non-permanent·es`))

2. Analyses de données spatiales : cartographie

Source : allison Horst @allison_horst

2.1. Notions de traitement de données spatiales : les objets, projections et données

Une organisation par couches

Les données sont organisées sous forme de couches superposables. On distingue généralement deux types de données : les données de type vecteur et les raster.

Les données vecteurs se définissent uniquement par des coordonnées. On trouvera des données vecteurs de type points, lignes et polygones. Un point sera défini par un couple de coordonnées XY, une ligne ou un polygone par les coordonnées de leurs sommets. Une couche vecteur sera soit de type point, soit de type ligne, soit de type polygone. Des données attributaires - supplémentaires - sont associées aux données spatiales : le nombre d’habitants d’une ville, par exemple.

Les données raster sont des images (vue satellite, carte administrative, carte IGN).

Projections et systèmes de coordonnées

La Terre est (à peu près) ronde et se présente en 3 dimensions. Afin de représenter les informations géographiques sur une carte en 2 dimensions, il faudra projeter les coordonnées spatiales, ce qui donnera lieu à des déformations. Les multiples méthodes de projections ont leurs avantages et inconvénients et seront utilisées selon différentes finalités.

Par exemple, la projection de Mercator conserver les angles mais pas les distances. Les surfaces et distances s’éloignant de l’équateur paraissent plus importantes qu’elles ne le sont en réalité. En revanche, cette projection conserve les angles et se trouve donc utile en navigation.

Extraits de “Le scandale de Mercator”, les Savoirs ambulants, https://www.lessavoirsambulants.fr/p/le-scandale-de-mercator.html

Un type de projection nous intéressera, il s’agit des projections coniques qui minimisent les déformations sur des zones du globes.

Utilisée en cartographie, on représente différentes zones en modifiant le référentiel. Le référentiel Lambert-93 est adapté pour représenter la France métropolitaine. Les systèmes de coordonnées sont multiples. Vous en rencontrerez deux principaux : WGS84 et Lambert-93 identifiables par les EPSG (European Petroleum Survey Group) 4326 et 2154.

Il est essentiel de s’assurer que tous vos fichiers spatiaux appartiennent au même système de coordonnées de référence (SCR). Les fichiers spatiaux peuvent être re-projetés automatiquement.

2.2. Les données spatiales dans R

sf est un des packages de référence pour traiter simplement des données spatiales dans R. Il a l’avantage d’être compatible avec l’esprit et les opérateurs tidyverse.

/!\ Comme tout le reste dans R, les choses évoluent, et rapidement. Notamment pour les packages les plus utilisés et maintenus ainsi que les formats de données les plus répandus.

Dès maintenant, installez les packages suivants (qui sont assez longs à installer, ne laissez pas cette ligne directement dans le script) :

install.packages(c("sf", "geojsonio","ggspatial"))

Puis, téléchargez le fichier “carto.zip” du moodle et mettez-le dans le dossier data :

Importer des données vecteurs

Exemple avec des premiers jeux de données sur la page d’exploration des données open data de la MEL :

Créer un dossier et extraire le contenu d’une archive zip :

# Créer un dossier "carto" dans le dossier data s'il n'existe pas déjà :
if(dir.exists("data/carto")==F){
dir.create("data/carto") 
}
# Extraire les fichiers du zip dans le dossier créé :
unzip(zipfile="data/carto.zip", exdir = "data/carto")

Le format shapefile

Shapefile = Le format shapefille est une des formes de stockage de données vectorielles qui s’est imposée (utilisée dans le SIG ArcGIS et compatible avec des logiciels libres comme QGIS). Il est en fait composé de plusieurs fichiers :

  • Un fichier .shp : Les données géographiques
  • Un .dbf : Base de données attributaires (informations sur les entités)
  • .shx : stockage des index
  • Eventuellement .sbn, .sbx (index de jointures), .xml (métadonnées).

Importer les fichiers spatiaux .shp

library(sf)
# Import du fichier shp :
garesMEL_sf <- st_read(dsn="data/carto/liste-des-gares.shp")
## Reading layer `liste-des-gares' from data source 
##   `/home/cecile/Documents/Cours et ateliers/ENSP_22-23/M2 - Approfondissement R/data/carto/liste-des-gares.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 59 features and 15 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: 2.804054 ymin: 50.52747 xmax: 3.234699 ymax: 50.76164
## Geodetic CRS:  WGS 84
class(garesMEL_sf)
## [1] "sf"         "data.frame"

Importer les fichiers spatiaux geojson

Vous trouverez aussi parfois des fichiers vectoriels en format GeoJSON. Pour les lire :

library(geojsonio)
# Importer un fichier GeoJSON :
communesMEL_sf <- geojson_read("data/carto/mel_communes.geojson", what="sp") %>%
  st_as_sf() # Le transformer en format sf

Représenter des données de vecteurs

Les objets spatiaux sf sont constitués d’une partie “géométrique” (“geometry”, informations spatiales) et d’une partie “données” (attributs qui seront représentés) sous la forme d’un data.frame.

names(communesMEL_sf)
## [1] "ut"           "nom"          "insee"        "territoire"   "surface"     
## [6] "code_posta"   "perimetre"    "geo_point_2d" "geometry"
names(garesMEL_sf)
##  [1] "code_uic"   "libelle"    "fret"       "voyageurs"  "code_ligne"
##  [6] "rg_troncon" "pk"         "commune"    "departemen" "idreseau"  
## [11] "idgaia"     "x_l93"      "y_l93"      "x_wgs84"    "y_wgs84"   
## [16] "geometry"

La partie geometry peut être représentée comme un type de graphique :

# L'objet `data` est précisé dans chaque geom_sf ou couche cartographique
ggplot()+
  geom_sf(data=communesMEL_sf)

Vérifier et modifier les systèmes de coordonnées

Il est important de vérifier le sytèmes de coordonnées : les graphiques reprojettent les couches à la volée mais certaines opérations ne seront pas possibles s’ils ne sont pas compatibles.

# st_crs(communesMEL_sf)
# st_crs(garesMEL_sf)

communesMEL_sf <- st_transform(communesMEL_sf,crs=2154)
garesMEL_sf <- st_transform(garesMEL_sf,crs=2154)

# st_crs(communesMEL_sf)
# st_crs(garesMEL_sf)

Représenter plusieurs couches sur une carte :

ggplot()+
  geom_sf(data=communesMEL_sf)+
  geom_sf(data=garesMEL_sf)

Représenter des données spatiales avec des informations à l’aide de ggplot2

La carte se construit donc par ajout d’éléments esthétiques, comme pour une graphique de données non-géographiques :

ggplot()+
  geom_sf(data=communesMEL_sf, aes(fill=surface))+
  viridis::scale_fill_viridis()+
  geom_sf(data=garesMEL_sf,aes(col=voyageurs))+
  scale_color_manual(values = c("black","grey"))+
  theme_bw()+
  labs(title = "Surface des communes et gares de la MEL",  col="Voyageurs",
       fill="Surface (m²)", caption="Source : MEL")+
  coord_sf( datum = NA)

Ajout d’attributs

Les données directement incluses dans les fichiers spatiaux ne sont pas toujours suffisantes et gagnet à être enrichies. Les éléments (points, polygones, lignes) possèdent souvent des identifiants. Par exemple, pour le fichier des communes de la MEL, la colonne “insee” correspond aux codes INSEE unique des communes (attentions, il ne s’agit pas des codes postaux).

Importer les données de la population légale des communes du recensement 2018 (fichier Communes.csv) et ajouter les informations au fichier spatial à l’aide d’une jointure de table avec left_join.

pop18 <- read.csv2("data/carto/Communes.csv") # import
communesMEL_sf <- left_join(communesMEL_sf, pop18, by =c("insee" = "DEPCOM")) # jointure
communesMEL_sf <- communesMEL_sf %>% filter(is.na(PTOT)==FALSE)

Représenter la population de chaque commune.

ggplot()+
  geom_sf(data=communesMEL_sf,aes(fill=PTOT),color="snow4",size=.2)+
  viridis::scale_fill_viridis()+
  theme_bw()+
  labs(title = "Population des communes de la MEL",  
       fill="Population", caption="Source : RP 2017, Insee")+
  coord_sf( datum = NA)

Représenter les communes selon les quintiles de population. On crée une variable représentant les quintiles de population avec les fonctions cut et quantile

communesMEL_sf$PTOT_quintiles <- cut(communesMEL_sf$PTOT,
                                     quantile(communesMEL_sf$PTOT, seq(0,1,.2), na.rm=T))

ggplot()+
  geom_sf(data=communesMEL_sf,aes(fill=PTOT_quintiles),
          color="snow4",size=.2)+
  viridis::scale_fill_viridis(discrete=T)+
  theme_bw()+
  labs(title = "Population des communes de la MEL",  
       fill="Population", caption="Source : RP 2017, Insee")+
  coord_sf( datum = NA)

Transformations des données spatiales

Au-delà de données contenues dans les fichiers spatiaux et de données extérieures les enrichissant, toute une série de transformation et de calculs peuvent être directement effectués sur les données spatiales : Combien y’a-t-il de gares dans chaque commune de la MEL ? Quel est le prix moyen d’achat des logements dans chaque quartier de Lille ?

Les fonctions de transformation et de calculs sur les fichiers spatiaux apparaissent avec le préfixe “st_”.

Importer les données DVF 2019

Les données DVF sont ” publiées et produites par la direction générale des finances publiques, permettent de connaître les transactions immobilières intervenues au cours des cinq dernières années sur le territoire métropolitain et les DOM-TOM, à l’exception de l’Alsace-Moselle et de Mayotte. Les données contenues sont issues des actes notariés et des informations cadastrales.” (https://www.data.gouv.fr/fr/datasets/5c4ae55a634f4117716d5656/)

Importer les données DVF (59.csv)

dvf_brut <- read.csv("data/carto/59.csv", stringsAsFactors = F)

Les données DVF sont des données issues des actes notariaux et pas forcément utilisables directement. Il faut passer par des traitements de nettoyage. Pour gagner du temps, j’ai préparé les données et fait les modifications suivantes :

  • En gardant uniquement les “nature_mutation” de type vente et vente en l’état futur d’achèvement ;
  • En gardant les ventes ayant des coordonnées géographiques et une valeur foncière renseignés ;
  • En imputant les valeurs manquantes des surface_reelle_bati et nombre_pieces_principales à 0 ;
  • En regroupant les lignes d’une même mutation (la valeur foncière est la même pour les différentes lignes, la surface et le nombre de pièces sont additionnés) ;
  • En gardant les 98 centiles de prix au m² centraux et les 98 centiles de surface bâties centraux ;
  • En ne gardant que les colonnes id_mutation, valeur_fonciere, surface_reelle_bati, nombre_pieces_principales, latitude, longitude

Le fichier nettoyé à importer est “dvf19_cours.rds”.

dvf <- readRDS("data/carto/dvf19_cours.rds")

Pour autant, ce n’est pas un fichier spatial … Mais on dispose de 2 colonnes contenant les informations “latitude” et “longitude”. On peut s’en servir pour créer un fichier spatial. On construit la partie “géométrique” à l’aide des colonnes longitude et latitude et de la fonction st_as_sf(). Les coordonnées sont exprimées selon la projection WGS84 (crs = 4326)

On peut transformer ce tableau de données en un fichier spatial à l’aide de st_as_sf en indiquant les colonnes de longitude (x) et latitude (y) et le crs (ici, 4326) :

dvf_sf <- dvf %>%
  st_as_sf(coords = c("longitude", "latitude"), crs=4326)
dvf_sf <- dvf_sf %>% st_transform(crs=2154)

Puis, on représente les points des ventes et les délimitations des communes de la MEL :

ggplot() +
  geom_sf(data=dvf_sf) +
  geom_sf(data=communesMEL_sf)

Pour la suite, on s’intéresse aux ventes immobilières à une échelle plus fine, celle des IRIS (unité territoriale infracommunale pour les communes de plus de 10.000 habitants) de la MEL (Métropole Européenne de Lille, fichier MEL_iris.shp). Importer le fichier spatial “MEL_iris” puis représentez les ventes immobilières à l’échelle des IRIS sur une carte.

irisMEL_sf <- st_read("data/carto/MEL_iris.shp")
## Reading layer `MEL_iris' from data source 
##   `/home/cecile/Documents/Cours et ateliers/ENSP_22-23/M2 - Approfondissement R/data/carto/MEL_iris.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 1271 features and 18 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 685041 ymin: 7044709 xmax: 719291.4 ymax: 7077561
## Projected CRS: RGF93_Lambert_93
# st_crs(MEL_iris)
ggplot() +
  geom_sf(data=dvf_sf) +
  geom_sf(data=irisMEL_sf)

Intersection

Pour associer chaque point à chaque iris, la fonction st_intersection permet de sortir la zone commune à deux fichiers sf en ajouter les données des iris au fichier des points des ventes immobilières.

dvf_MEL_sf <- st_intersection( dvf_sf, irisMEL_sf)
ggplot() +
  geom_sf(data=dvf_MEL_sf)+
  geom_sf(data=irisMEL_sf, fill="transparent")

On a oublié de regarder les systèmes de projection de ces 2 objets ! Ayez ce réflexe, même si les fonctions de carto reprojetent les données à la volée.

Compter les transactions dans chaque IRIS

Un fichier sf se compose d’un tableau de données et des informations spatiales associées à chaque entité.

Produire une table contenant le nombre de transactions immobiliers par IRIS (CODE_IR) à l’aide des fonctions group_by et summarise.

transactions_IRIS <- dvf_MEL_sf %>% 
  group_by(CODE_IR) %>% 
  summarise(nb_transactions=n()) %>% 
  data.frame() %>% 
  select(CODE_IR, nb_transactions)

Représenter cette données sur une carte (le fonds de carte des IRIS est irisMEL_sf)

irisMEL_sf <- left_join(irisMEL_sf, transactions_IRIS, by="CODE_IR")

ggplot()+ geom_sf(data=irisMEL_sf,size=.2,
                  aes(fill=nb_transactions))+
  viridis::scale_fill_viridis()+
  theme_bw()+
coord_sf(crs = st_crs(irisMEL_sf), datum = NA)

Exercice

Représenter la taille moyenne (variable surface_reelle_bati) des logements vendus par IRIS dans la MEL puis à Lille.

taille <- dvf_MEL_sf %>% 
  # filter(surface_reelle_bati<800) %>% 
  group_by(CODE_IR) %>% 
  summarise(taille = mean(surface_reelle_bati)) %>% 
  data.frame() %>%
  select(CODE_IR, taille)

irisMEL_sf$taille <- NULL
irisMEL_sf <- left_join(irisMEL_sf, taille, by="CODE_IR")

ggplot()+
  geom_sf(data=irisMEL_sf, size=.1)+
  aes(fill=taille)+
  viridis::scale_fill_viridis()+
  theme_bw()+
  labs(title = "Taille des logements vendus dans la métropole lilloise",
       fill="Surface réelle du bâti (m²)",
       caption="Source : DVF, Ministère des Finances, 2019")

irisLille_sf <- irisMEL_sf %>% filter(NOM_COM=="Lille")

ggplot()+
  geom_sf(data=irisLille_sf, size=.1)+
  aes(fill=taille)+
  viridis::scale_fill_viridis()+
  theme_bw()+
  labs(title = "Taille des logements vendus à Lille",
       fill="Surface réelle du bâti (m²)",
       caption="Source : DVF, Ministère des Finances, 2019")

Représenter les lieux des ventes dans Lille et les prix au m² de ces ventes.

dvf_Lille_sf <-  st_intersection( dvf_MEL_sf, irisLille_sf)
dvf_Lille_sf$prix <- round(dvf_Lille_sf$valeur_fonciere/dvf_Lille_sf$surface_reelle_bati)

ggplot()+
  geom_sf(data=irisLille_sf, size=.1, fill="transparent")+ 
  geom_sf(data=dvf_Lille_sf, aes(col=prix), size=.1)+
  viridis::scale_color_viridis()+
  coord_sf(crs = st_crs(irisLille_sf), datum = NA)+
  theme_bw()

2.3. Ajout de raster ou d’éléments de fond de carte

Raster

Le package ggspatial permet d’ajouter un fond de type raster (image ou “tuiles”).

library(ggspatial)
 
 ggplot()+
  annotation_map_tile(zoom=13,type="cartolight")+
  geom_sf(data=irisLille_sf, size=.3, fill="transparent", col="grey60")+
  geom_sf(data=dvf_Lille_sf, aes(col=prix), size=.1)+
  viridis::scale_color_viridis()+
  theme_bw()+
  labs(title = "Prix au m² des ventes de 2019 à Lille",
       col="Prix (/m²)",
       caption="Sources : DVF, Ministère des Finances, 2019, fond : OSM")+
coord_sf(crs = st_crs(irisLille_sf), datum = NA)

Extraire des éléments de fond d’OpenStreetMap

Si on souhaite faire apparaître certains éléments de fond au lieu des tuiles images, on peut faire une extraction depuis OpenStreetMap à l’aide du package {osmdata} :

library(osmdata)
roads <- opq(bbox = st_bbox(irisLille_sf) ) %>%
  add_osm_feature(key = 'highway') %>%
  osmdata_sf()
roads <- roads$osm_lines
st_crs(roads) <- 4326 
roadsLille_sf <-  st_intersection( roads, irisLille_sf)

ggplot()+geom_sf(data=roads)
ggplot()+geom_sf(data=roadsLille_sf)

Cartes interactives

Enfin, les cartes peuvent être rendues interactives ; C’est notamment utile si vous en fait une application {shiny} (voir séminaire/TD dédié). Plusieurs packages le permettent, vous pouvez vous reporter à ces deux articles :

Un exemple avec {tmap} :

library(tmap)
library(viridis)
tmap_mode(mode = "view") # Active le mode interactif
tm_basemap("CartoDB.Positron")+ # tuiles rasters
  tm_shape(dvf_Lille_sf)+ # données
  tm_symbols(col = "prix", border.col = "transparent", 
             scale = .3, alpha = 1, palette =viridis(n=10) ,
             perceptual = TRUE) # Taille des points reste la même

Un dernier exercice

Représenter de la façon la plus lisible le prix au m² selon les catégories de taille des ventes de la MEL. Essayez de trouver une représentation graphique originale et adaptée à ce que vous voulez raconter !

# Prix au m² :
dvf_MEL_sf$prix <- round(dvf_MEL_sf$valeur_fonciere/dvf_MEL_sf$surface_reelle_bati)

# Catégories de tailles :
dvf_MEL_sf$taille <- case_when(dvf_MEL_sf$surface_reelle_bati<40 ~ "40-",
                               dvf_MEL_sf$surface_reelle_bati<=60 ~ "40-60",
                               dvf_MEL_sf$surface_reelle_bati<=80 ~ "60-80",
                               dvf_MEL_sf$surface_reelle_bati<=100 ~ "80-100",
                               dvf_MEL_sf$surface_reelle_bati<=120 ~ "100-120",
                               dvf_MEL_sf$surface_reelle_bati>120 ~ "120+",
                               T~"")   
dvf_MEL_sf$taille <- factor(dvf_MEL_sf$taille,
                              levels=c("40-","40-60","60-80","80-100",
                                       "100-120","120+"))

# Graphique 1 :

# TaillePrix
 TaillePrix <- dvf_MEL_sf %>% 
  group_by(taille) %>% 
  summarise("Prix moyen"=mean(prix),
            "Prix médian"=median(prix))

TaillePrix %>%  
  ggplot()+ 
  geom_linerange(aes(ymax = `Prix moyen`,
                     ymin=`Prix médian`,
                     x=taille), col="grey60", size=1.5)+
  geom_point(aes(y= `Prix moyen`, x=taille), col="white",
            size=1.5)+
  geom_point(aes(y= `Prix moyen`, x=taille), col="firebrick",
             shape=1, stroke=2, size=1.5)+
  geom_point(aes(y= `Prix médian`, x=taille), col="white",
             size=1.5)+
  geom_point(aes(y= `Prix médian`, x=taille), col="seagreen",
             shape=1, stroke=2, size=1.5)+
  theme_bw()+ coord_flip()+ 
  labs(y="Prix (rouge=moyen, vert=médian)",
       x="Taille (m²)",
       title="Prix/m² des logements selon la taille",
       subtitle="Métropole Européenne de Lille (2019)",
       caption="Source : DVF")

# Graphique 2 :
library(ggridges)
dvf_MEL_sf %>% ggplot(aes(x=prix, y=taille, fill=stat(x)))+
  geom_density_ridges_gradient(col="white",scale = 3, size = .5, rel_min_height = 0.01) +
  scale_fill_viridis_c(name = "", option = "C")+
  theme_bw()

# https://www.datanovia.com/en/blog/elegant-visualization-of-density-distribution-in-r-using-ridgeline/

3. Données pondérées

Retour aux données Coco sur lesquelles vous allez travailler pour le rendu de fin de semestre :

Il se trouve que les données Coco sont des données construites sur échantillon. On l’a déjà entre-aperçu, la structure de la population est particulière alors que le panel se veut représentatif de la population de France métropolitaine. Pour corriger les effets de structure de l’échantillon enquêté, les producteur·rices des données ont construit des coefficients de redressement que nous allons utiliser.

On a vu que le panel était plutôt âgé. On retrouve cette correction dans les poids :

coco %>% 
  group_by(age_quinq) %>% 
  summarise(`poids médian`=median(POIDS),
            `poids moyen`=mean(POIDS))

Les plus jeunes sont “sur-pondéré·es” et les plus âgé·es “sous-pondéré·es”. Et la somme des POIDS est égale à la taille de l’échantillon, on dit que la pondération est normée.

Plusieurs packages permettent d’utiliser les pondérations :

  1. Pour les traitements de types “tris à plat” et “tris croisés”, on utilisera la fonction très pratique wtd.table du package {questionr}.
  2. {lourdR} : un package de Cécile Rodrigues (@Grisoudre), pas le plus convivial et un choix de nom discutable mais pas si inutile que ça…
  3. Pour d’autres traitements numériques, nous verrons {survey} qui est en général présenté dans les tutoriels et enseignements de R en SHS (https://larmarange.github.io/analyse-R/donnees-ponderees.html). {survey} est plus délicat à prendre en main que d’autres packages mais il intègre l’ensemble des traitements et tests qu’on peut vouloir faire. Il faut parfois un peu “bidouiller” des tableaux de résultats pour les présenter.

3.1. Tris à plat et tris croisés avec {questionr}

On commence par faire le tri à plat redressé des PCS. Il faut déjà recoder les PCS de la variable coco1_q35 :

coco$pcs <- fct_collapse(factor(coco$q35),
                       "Agriculteurs"="1",
                       "Artisans, comm., chef d'E"="2",
                       "Cadres et pr.intell.sup."="3",
                       "Prof. intermédiaires"="5",
                       "Employés"="4","Ouvriers"="6",
                       NULL = c("7","9999"))
freq(coco$pcs)

Effectifs

# N pondéré :
wtd.table(coco$pcs, weights = coco$POIDS,
          useNA = "ifany")
##              Agriculteurs Artisans, comm., chef d'E  Cadres et pr.intell.sup. 
##                  26.78782                  82.68868                 222.98026 
##                  Employés      Prof. intermédiaires                  Ouvriers 
##                 167.20659                 448.01096                  87.77051 
##                      <NA> 
##                  40.55518

Du fait des poids, les effectifs ont des décimales donc on les arrondit :

wtd.table(coco$pcs, weights = coco$POIDS,
          useNA = "ifany") %>% 
  round()#<<
##              Agriculteurs Artisans, comm., chef d'E  Cadres et pr.intell.sup. 
##                        27                        83                       223 
##                  Employés      Prof. intermédiaires                  Ouvriers 
##                       167                       448                        88 
##                      <NA> 
##                        41

Pourcentages

# % pondérés :
round(wtd.table(coco$pcs, weights = coco$POIDS,
          useNA = "ifany") %>% 
  prop.table()*100,1)#<<
##              Agriculteurs Artisans, comm., chef d'E  Cadres et pr.intell.sup. 
##                       2.5                       7.7                      20.7 
##                  Employés      Prof. intermédiaires                  Ouvriers 
##                      15.5                      41.6                       8.2 
##                      <NA> 
##                       3.8

Tris croisés

Croisons la PCS avec la variable q38 (Situation actuelle vis-à-vis du télétravail) :

# Recodages :
coco$q38 <- fct_recode(factor(coco$q38),
                               "Lieu habituel"="1",
                               "Alternance lieu habituel / TT"="2",
                               "Télé-travail"="3", NULL="6666")
# Tableau :
wtd.table(coco$pcs, coco$q38, 
          weights = coco$POIDS, useNA = "ifany") %>% 
  rprop(digits=1)
##                           Lieu habituel Alternance lieu habituel / TT
## Agriculteurs               42.2           3.0                        
## Artisans, comm., chef d'E   2.8           1.0                        
## Cadres et pr.intell.sup.    4.4           3.3                        
## Employés                   14.0           2.6                        
## Prof. intermédiaires       18.5           1.5                        
## Ouvriers                   12.0           0.0                        
## <NA>                        0.0           0.0                        
## Ensemble                   13.1           1.8                        
##                           Télé-travail <NA>  Total
## Agriculteurs                0.0         54.7 100.0
## Artisans, comm., chef d'E   4.5         91.6 100.0
## Cadres et pr.intell.sup.   27.6         64.6 100.0
## Employés                   12.2         71.2 100.0
## Prof. intermédiaires        5.6         74.4 100.0
## Ouvriers                    0.8         87.2 100.0
## <NA>                       10.9         89.1 100.0
## Ensemble                   10.8         74.3 100.0

3.2. Résultats pondérés avec {lourdR}

{lourdR} est un package que j’ai commencé à faire. Il est perfectible mais utile. Je l’ai déposé sur github, pour l’installer, voici les lignes à exécuter (une fois) :

if (!require("devtools")) install.packages("devtools", dep=T)
# installe devtools si vous ne l'avez pas
devtools::install_github("Grisoudre/lourdR") # installe lourdR

Ensuite, il peut être chargé comme un package classique :

library(lourdR)

{lourdR} contient - pour le moment - deux fonctions, on utilise déjà freqp : copie de la fonction freq présentant les résultats non-pondérés et pondérés :

freqp(coco, var = "pcs", poids = "POIDS")
# n=F : sans les effectifs
# p=F : sans les pourcentages
# brut=F : sans les valeurs brutes
# exclude = "item" : Valeur à exclure

On peut voir l’effet des redressements sur les résultats !

L’autre fonction freqm présente dans un même tableau des colonnes avec un même préfixe et des mêmes modalités - typique des variables issues de questions à choix multiples. Prenons les variables q24_1 à q24_7 (Vit avec : …) :

freqm(coco, "q24", calc="%")
# calc = "n", "%", ou "all"
# exclude = "item" : valeur exclue
# poids = "variable_poids" : si données pondérées

Puis on renomme les modalités de “name”. On peut également pondérer ce tableau.

En retirant les “non-concerné·es”, en pondérant :

freqm(coco, "q24", calc="%", exclude="6666", poids="POIDS")

Un peu mieux mis en forme avec {tidyverse} :

freqm(coco, "q24", calc="%", exclude="6666", poids="POIDS") %>% 
  mutate(name= fct_recode(name,
                          "Conjoint·e"="q24_1",
                          "Enfants/petits-enfants"="q24_2",
                          "Parents/Grands-parents"="q24_3",
                          "Famille (autre)"="q24_4",
                          "Ami·es ou autres"="q24_5",
                          "Chien(s)"="q24_6",
                          "Autre(s) animal(aux)"="q24_7")) %>% 
  rename("Vit avec"="name", "Oui"="X1") %>% 
  select(-X2, -Total) %>% 
  arrange(desc(Oui))

3.3. Résultats pondérés avec {survey}

Pour le moment, on a utilisé les poids pour présenter des résultats de types effectifs et pourcentages sur effectifs. C’est en fait assez limité. Si on souhaite faire des calculs plus avancés - types régressions - on peut se référer à un autre package, {survey}. Il est en général présenté dans les tutoriels et enseignements de R en SHS (https://larmarange.github.io/analyse-R/donnees-ponderees.html). Il est plus délicat à prendre en main que d’autres packages mais il intègre l’ensemble des traitements et tests qu’on peut vouloir faire. Il faut parfois un peu “bidouiller” des tableaux de résultats pour les présenter. Installons-le :

install.packages("survey", dep = T)
library(survey)

Avant toute chose, il faut créer un nouvel objet de type “survey design” explicitant le plan d’échantillonnage (ou design) et donc la variable de pondération à considérer. Donc on produit le survey design avec svydesign :

coco$sexe <- fct_recode(factor(coco$eayy_A1),"Homme"="1","Femme"="2")
# Recodage du genre avant
cocow <- svydesign(data = coco,
                        weights = ~ coco$POIDS,
                        ids = ~1)
class(cocow) # = Objet bizarre
## [1] "survey.design2" "survey.design"

L’objet cocow n’est plus un tableau de données. Le tableau de données se trouve en fait dans la partie “variables” du “survey.design” :

names(cocow$variables)[1:10]
##  [1] "UID_ANONYM_COVID" "q01"              "q02"              "q03"             
##  [5] "q04"              "q05_c"            "q06_1"            "q06_2"           
##  [9] "q06_3"            "q06_4"

Donc !!!!! Il est important de faire tous vos recodages et sélections d’individus AVANT de créer l’objet survey design. N’oubliez pas qu’un script est linéaire.

On peut toujours faire des calculs sur le tableau initial :

freq(cocow$variables$sexe)

On peut aussi être tenté·e de faire les recodages et modifications sur ce tableau-là (je le fais aussi) mais il faut être sûr·e de soi, ça complique vite les choses. Conseil : Faire les recodages, créations de variables, etc. avant !

Effectifs pondérés

  • Compter des effectifs pondérés : svytable(~variable, design)
svytable(~sexe, cocow)
## sexe
##    Homme    Femme 
## 529.2982 546.7018

–

  • En arrondissant les effectifs :
svytable(~sexe, cocow, round = T)
## sexe
## Homme Femme 
##   529   547
  • Présenté dans un tableau (éventuellement exportable via freq) :
svytable(~sexe, cocow, round = T) %>% 
  freq(valid=F, total = TRUE)

!! à la gestion des NA

svytable(~pcs, cocow, round=T) %>% 
  freq(valid=F, total = TRUE)

{survey} ignore les NA !!

–

Il faut les expliciter

cocow$variables$pcs <- fct_explicit_na(cocow$variables$pcs,"NR")
svytable(~pcs, cocow, round=T) %>% 
  freq(valid=F, total = TRUE)

Tableaux croisés pondérés

svytable(~pcs+q38, cocow) %>% 
  rprop(digits=1)
##                            q38
## pcs                         Lieu habituel Alternance lieu habituel / TT
##   Agriculteurs               93.3           6.7                        
##   Artisans, comm., chef d'E  33.7          12.2                        
##   Cadres et pr.intell.sup.   12.6           9.3                        
##   Employés                   48.6           8.9                        
##   Prof. intermédiaires       72.3           5.8                        
##   Ouvriers                   93.5           0.0                        
##   NR                          0.0           0.0                        
##   Ensemble                   50.8           7.2                        
##                            q38
## pcs                         Télé-travail Total
##   Agriculteurs                0.0        100.0
##   Artisans, comm., chef d'E  54.1        100.0
##   Cadres et pr.intell.sup.   78.2        100.0
##   Employés                   42.4        100.0
##   Prof. intermédiaires       21.9        100.0
##   Ouvriers                    6.5        100.0
##   NR                        100.0        100.0
##   Ensemble                   42.0        100.0

Autres calculs sur des variables numériques

Exemple d’une moyenne pondérée (q12 = nombre de personnes dans l’entourage ayant contracté le covid) :

svymean(~ q12, cocow, na.rm=T)
##       mean     SE
## q12 2.8041 0.1465

3.4. Graphiques de données pondérées

{ggplot2} permet d’intégrer les données de survey design. Il faut ajouter deux précision :

  • La table ‘variables’ du survey design dans la partie ggplot()
  • Le poids du survey design dans la partie aes()

C’est un peu du bidouillage…

Ex : Graphique du lieu d’activité professionnelle pendant le 1er confinement selon la PCS :

# Graphique :
ggplot(cocow$variables) +#<<
  aes(weight = weights(cocow),#<<
      x = pcs, fill = q38) +
  geom_bar(position = "fill")+
  theme(axis.text.x = element_text(angle=45, hjust=1))+
  labs(title = "Lieu de travail au cours des deux dernières semaines",
       x="PCS",fill="",y="",
       caption ="Source : Coco1 - ELIPSS, CDSP, avril 2020")

  • Exercices supplémentaires facultatifs : Refaire toutes les questions de la partie I. en version pondérée

  • Aller plus loin pour traiter les données pondérées avec {survey} : https://larmarange.github.io/analyse-R/donnees-ponderees.html

  • Nous verrons la façon de produire des tests statistiques sur des données pondérées dans la partie suivante.

4. Statistiques inférentielles

4.1. Eléments de vocabulaire

Echantillon :

= Une partie de la population qu’on souhaite interroger. Les personnes interrogées – l’échantillon – sont choisies de façon à être représentatives de l’ensemble de la population : La manière dont sont choisies ces personnes est appelée échantillonnage.

Représentativité :

La représentativité statistique permettra de généraliser des résultats et observations faites sur l’échantillon à l’ensemble de la population. Pour être représentatif, un échantillon doit présenter les mêmes caractéristiques que sa population (soit par règle statistique d’aléatoire, soit par construction “raisonnée” d’un échantillon).

La représentativité statistique s’appuie sur les notions de probabilité. Elle permettra d’adopter un raisonnement déductif : C’est à dire de partir du particulier (l’échantillon) pour généraliser des résultats observés à un ensemble plus vaste (la population) = Statistiques inférentielles.

La statistique inférentielle consiste à généraliser des résultats observés sur un échantillon à l’ensemble de la population de référence. Produire des statistiques inférentielles se fait donc sous certaines conditions. Une des premières conditions nécessaires à la réalisation de statistiques inférentielles porte sur l’échantillon ; celui-ci doit pouvoir garantir une certaine représentativité statistique :

  • Echantillons probabilistes : chaque individu de la population doit avoir une probabilité non-nulle et connue de faire partie de l’échantillon. Ces échantillons supposent d’avoir une base de sondage (une liste exhaustive des individus constituant la population-mère), ce qui est rare.
  • Echantillons empiriques : Autres types d’échantillons, plus ou moins raisonnés, comme les échantillons par quotas par exemple. L’échantillon par quotas a une structure similaire à celle de la population-mère selon certains critères (sexe, âge, CSP, etc.).

4.2. Estimation d’une moyenne

Avant de traiter les tests statistiques, voici deux exemples d’estimation de valeurs. Dans cette partie, à partir de données fictives, on découvrira comment estimer une valeur dans la population générale en utilisant qu’une partie des individus de cette population (un échantillon).

Exemples pour observer comment le choix de l’échantillon agit sur les résultats :

Soit, pour l’exemple 1, une population connue fictive de 10 000 individus qui ont mangé du chocolat l’année dernière.

Données fictives

On construit (fictivement) un tableau reprenant la quantité (en kg) de chocolat mangée l’an dernier par chacune de ces 10 000 personnes au cours de l’année écoulée. toutes ces données sont habituellement inconnues (c’est la quantité de chocolat consommée par toute la population d’un pays, par exemple, qu’on ne peut pas interroger de façon exhaustive).

dt <- data.frame(choco=round(rnorm(10000, mean = 13.2 , sd = 5)))
dt[dt<0]<-0
# loi normale de moyenne 13,2 et d'écart type 5

head(dt)

L’individu avec l’identifiant 1 a mangé 21 kg de chocolat, la personne avec l’identifiant 2 en a mangé 11 kg.

La moyenne (habituellement inconnue et qu’on veut estimer) de cette quantité de chocolat ingurgitée est connue et de :

mean(dt$choco)
## [1] 13.1801

= Valeur qu’on veut estimer !

Et les individus se répartissent comme ceci selon la quantité de chocolat mangée pendant l’année :

ggplot(dt)+
  aes(x=choco) + 
  geom_histogram(stat="count")

On reconnaît la forme de la répartition normale.

Construction d’un échantillon

On demande à 500 d’entre elles et eux, l’échantillon, choisi·es au hasard (aléatoire simple), la quantité de chocolat mangée au cours des 12 derniers mois. On établit ensuite la moyenne de chocolat consommé par ces 500 individus, soit la moyenne estimée.

  • Sélection aléatoire de 500 individus dans la population et calcul des statistiques de base (dispersion, moyenne, écart-type) :
dt$sample <- sample(1:10000,10000,replace=F)

tirage <- data.frame(tirage=1, 
                     moyenne=mean(dt[dt$sample <501,"choco"]),
                     sd = sd(dt[dt$sample <501,"choco"]),
                     n=nrow(dt[dt$sample <501,]))
tirage

Sur ce 1er tirage, on a une moyenne de 13.27.

On répète une 2ème fois l’opération avec 500 individus choisis au hasard parmi la population de départ ; on obtient une moyenne de :

dt$sample <- sample(1:10000,10000,replace=F)

tirage2 <- data.frame(tirage=2, 
                     moyenne=mean(dt[dt$sample <501,"choco"]),
                     sd = sd(dt[dt$sample <501,"choco"]),
                     n=nrow(dt[dt$sample <501,]))
tirage <- rbind(tirage, tirage2)
tirage

On a cette fois une moyenne de 13.312.

Répéter l’action avec une boucle

On répète l’opération encore 498 fois. On a donc tiré 500 échantillons constitués de 500 personnes. Pour éviter de taper ou exécuter le script 498 fois, on la produit dans une boucle : soit un i allant de 2 à 500, pour chacune des valeurs, on reproduit un tirage aléatoire nouveau et on ajoute la ligne i au tableau des tirages avec le numéro de l’échantillon, la moyenne calculée, l’écart-type calculé et le nombre d’individus dans l’échantillon.

for (i in 2:500) # "éléments" sur lesquels va tourner la boucle
  { # ouverture de la boucle
dt$sample <- sample(1:10000,10000,replace=F) # ajout d'un nombre aléatoire entre 1 et 10000 sans remise
tirage[i,] <- c(i,  # Ajout de la ligne i avec le nombre i
                mean(dt[dt$sample <501,"choco"]), # la moyenne du tirage de 500
                sd = sd(dt[dt$sample <501,"choco"]), # l'écart-type du tirage de 500
                n=nrow(dt[dt$sample <501,])) # le nombre de lignes dans le tirage
}
head(tirage)

Histogramme de l’ensemble des moyennes de ces 500 échantillons :

ggplot(tirage)+
  aes(x=moyenne, stat="count") + 
  geom_histogram( bins = 20)

Sur la base d’un tirage aléatoire simple, en fonction de l’échantillon choisi, la moyenne de la quantité de chocolat mangé l’année dernière varie (ici, assez peu). Cette variabilité est due au hasard de l’échantillon ; on observe une marge d’erreur.

Les conclusions sur la population comportent toujours une marge d’incertitude. Pour être sûr·es à 100% qu’un résultat soit vrai, il faut un échantillon exhaustif ; ie interroger toute la population des 10 000 de départ. Ce n’est que rarement possible. On cherche à probabiliser l’incertain.

On considère, en s’appuyant sur la loi des grands nombres, qu’on s’approche des caractéristiques de la population en augmentant la taille de l’échantillon (tout en ayant l’hypothèse nécessaire de l’aléatoire de l’échantillon). On calcule une moyenne, une proportion et on peut en déduire la probabilité qu’un intervalle contiennent la moyenne ou la proportion pour l’ensemble de la population :

Cette probabilité correspond au seuil de confiance qu’on décide de se donner. L’intervalle contenant – probablement – la moyenne pour la population est l’intervalle de confiance. La valeur calculée sur l’échantillon est finalement une estimation. Note : Il n’est pas sûr que la vraie valeur se trouve dans l’intervalle, ce n’est pas une certitude absolue.

Cette estimation sera en fait un intervalle : On dira qu’« On estime qu’il y a une probabilité de 95% que la moyenne de chocolat consommé chaque année par la population française se situe dans l’intervalle [12,7 ; 13,6].»

Intervalle de confiance

L’intervalle de confiance de la moyenne calculée peut être estimé grâce à l’écart-type (dans l’échantillon, la dispersion de la quantité de chocolat consommée l’an dernier est-elle élevée ou non ?). L’estimation à 95% de la moyenne sera l’intervalle :

[ x̅ - (1,96 σ / √n) ; x̅ + (1,96 σ / √n ) ]

avec x̅ : Moyenne ; σ : écart-type ; n : taille de l’échantillon.

Pour le 1er échantillon tiré : Si notre écart-type est de 5.16 l’estimation de la moyenne sera : [ 13.15 – (1.96 x 5.16 / √500) ; 13.15 + (1.96 x 5.16 / √500) ] = [13.15 - ( 1.96 x 5.16 / 22.36) ; 13.15 + ( 1.96 x 5.16 / 22.36) ] = [ 13.15 – 0.45 ; 13.15 + 0.45] = [ 12.7 ; 13.6 ]

On peut dire qu’on a 95 % de chance que la quantité de chocolat consommé soit situé entre :

tirage$moyenne[i] - (1.96 * tirage$sd[i] / (tirage$n[i]^.5) )
## [1] 12.34567

et

tirage$moyenne[i] + (1.96 * tirage$sd[i] / (tirage$n[i]^.5) )
## [1] 13.20233

La marge d’erreur étant de :

(1.96 * tirage$sd[i] / (tirage$n[i]^.5) )
## [1] 0.4283278

et la moyenne calculée de :

mean(tirage$moyenne[i])
## [1] 12.774

Même procédure avec un échantillon de taille 50 individus, cette fois-ci

i<-501
dt$sample <- sample(1:10000,10000,replace=F)
tirage[i,] <- c(i, 
                mean(dt[dt$sample <51,"choco"]),
                sd = sd(dt[dt$sample <51,"choco"]),
                n=nrow(dt[dt$sample <51,]))

tirage$moyenne[i] - (1.96 * tirage$sd[i] / (tirage$n[i]^.5) )
## [1] 11.6418
tirage$moyenne[i] + (1.96 * tirage$sd[i] / (tirage$n[i]^.5) )
## [1] 14.2382
(1.96 * tirage$sd[i] / (tirage$n[i]^.5) )
## [1] 1.298198
mean(dt$choco)
## [1] 13.1801

L’intervalle de confiance (et donc la précision de l’estimation) est influencé par :

  • Le niveau de confiance -> au choix, arbitraire, se reporter ensuite à la table de la loi correspondante (ici, loi normale)

  • La taille de l’échantillon (racine de n) -> suivant les moyens de l’enquête ; plus l’échantillon est grand, plus l’intervalle de confiance est petit

  • La dispersion des valeurs (écart-type) -> indépendant du chercheur, nécessité de l’observer. Plus il est faible, plus l’intervalle est petit.

Différence entre des échantillons aléatoires constitués de 50 et de 500 individus :

for (i in 502:1000){
dt$sample <- sample(1:10000,10000,replace=F)
tirage[i,] <- c(i, 
                mean(dt[dt$sample <51,"choco"]),
                sd = sd(dt[dt$sample <51,"choco"]),
                n=nrow(dt[dt$sample <51,]))
}


ggplot(tirage)+aes(x=moyenne,stat="count") + 
  geom_histogram() +
  stat_bin(bins=30)+ 
  facet_grid(~n)

Donc, sur un échantillon donné, les calculs qui seront produits seront assortis d’une marge d’erreur. Il s’agit d’une estimation car l’échantillon reste un échantillon parmi d’autres possibles.

4.3. Estimation d’une proportion

En sciences sociales, on traite surtout des données nominales. Et donc on recherche des proportions (quelle proportion des enquêté·es aiment le chocolat). Prenons donc un deuxième exemple : L’institut de sondage OPIF s’intéresse au plat du dimanche préféré des habitant·es d’un pays. Après un 1er tour élargi, deux plats-candidats restent en lice pour le 2nd : le Dahl de lentilles et le Boeuf bourguignon.

Données fictives

dt2 <- data.frame(repas_pref= c(rep("Dahl de lentille",5100), rep("Boeuf bourguignon", 4900)))

La proportion réelle de vote pour le Dahl de lentilles est de 51%. On veut l’estimer en prenant un échantillon de cet ensemble. En reproduisant la démarche du 1er exemple, on tire 1.000 échantillons différents pour voir les variations : le Dahl de lentilles va-t-il l’emporter ?

tirage <- NULL
for (i in 1:1000){
  dt2$sample <- sample(1:10000,10000,replace=F)
 tirage <- rbind(tirage,
                  data.frame(tirage=i, 
                     prop=round(nrow(dt2[dt2$sample<501 & dt2$repas_pref=="Dahl de lentille",])/
                                  nrow(dt2[dt2$sample<501,]),3)))
}

Graphique des estimations pour ces différents tirages :

ggplot(tirage)+aes(x=prop)+ geom_histogram()+geom_vline(xintercept=.5, col="red")

A priori, plus de tirages semblent indiquer que le Dahl l’emporterait mais pas tous. En estimant les résultats du Dahl de lentilles avec un seuil de confiance de 95%, on obtient l’intervalle de confiance suivant :

Intervalle de confiance

L’intervalle de confiance d’une proportion s’appuie sur la proportion calcullée et la taille de l’échantillon. L’estimation à 95% de la proportion sera l’intervalle :

[ P - 1.96 * √(P (1 - P) / n) ; P + 1.96 * √(P (1 - P) / n)]

avec P : Probabilité calculée ; n : taille de l’échantillon.

Pour le 1er échantillon tiré : Si notre proportion calculée est de 0.52 l’estimation de la proportion sera :

p=.52
ic=  1.96*(((p*(1-p))/500)^.5)
print(paste0("[",round(p-ic,3)," - ",round(p+ic,3),"]"))
## [1] "[0.476 - 0.564]"

Si on souhaite estimer cette proportion avec une probabilité de 90% et qu’on interroge 2.000 personnes :

p=.52
ic=  1.64*(((p*(1-p))/2000)^.5)
print(paste0("[",round(p-ic,3)," - ",round(p+ic,3),"]"))
## [1] "[0.502 - 0.538]"

On pourrait alors dire que le Dahl a toutes les chances de l’emporter (avec une probabilité de 90%).

4.4. Tests statistiques

Les tests statistiques permettent d’établir si les résultats obtenus sur l’échantillon sont dus au hasard ou s’ils se rapprochent du résultat qu’on aurait observé sur la population. Ils sonts mobilisés pour confirmer la probabilité (ou plausibilité) d’une corrélation ou d’un lien entre deux variables. La plupart du temps, ils permettent de calculer le risque de se tromper en établissant une hypothèse et ont pour objectif de connaître le niveau de signification (fiabilité) d’un résultat.

On produit une hypothèse appelée l’hypothèse nulle (H0) qui va supposer que les résultats obtenus sont dus au hasard et donc les variables observées sont indépendantes. Le test vise à savoir si on peut rejeter H0 pour affirmer que H1 (hypothèse inverse à H0) est probablement vraie et que les deux variables sont liées.

Pour réaliser et valide ou non ce type d’hypothèses, on va :

  • Faire une analyse descriptive : observer les résultats sur l’échantillon ;
  • Tester la significativité du résultat obtenu : soit la probabilité de se tromper en estimant que le résultat observé est valable pour l’ensemble de la population mère ;
  • En fonction du résultat du test, on pourra établir une interprétation sociologique de la corrélation observée.

Test d’indépendance du khi2

A partir d’un tableau d’effectifs croisés, si les proportions semblent montrer que les 2 variables sont liées (jouent dans le même sens ou inversement), le test d’indépendance du khi-2 permet de tester s’il est probable ou non que cette observation soit due au hasard.

Le khi-deux, chi-carré ou chi-square, mesure les écarts entre des effectifs observés et les effectifs théoriques (distribution aléatoire) à partir des marges du tableau croisé (appelé tableau de contingence). Le test du khi2 répond à la question : « Quelle est la probabilité pour que les deux variables qualitatives soient indépendantes l’une de l’autre ? Soit pour que leur distribution se fasse au hasard et qu’il n’y ait aucune relation statistique entre elles ». Il permet d’établir à quel seuil de probabilité on peut rejeter l’hypothèse d’indépendance entre des effectifs observés et les effectifs théoriques. Il nous donne une information sur la corrélation entre deux variables qualitatives.

Un test de khi deux s’applique uniquement sur des tableaux de contingence :

  • ayant au moins 2 lignes et 2 colonnes,
  • contenant des valeurs positives entières,
  • ayant au moins 60 observations au total,
  • ayant au minimum 5 observations par cases du tableau et/ou dans le tableau des effectifs théoriques.

Exemple : A partir des données Coco1

On souhaite savoir si le changement de situation professionnelle dû à la pandémie peut être relié à la pcs des enquêté·es. D’après les données de l’enquête, est-ce qu’on peut dire que la catégorie socio-professionnelle a un effet sur le fait d’avoir poursuivi son activité, été mis en congé ou au chômage ? Nous avons déjà créé la table coco.activite conservant la sous-population des enquêté·es en activité avant la crise sanitaire ainsi que la variable sitpro.conf de la situation professionnelle pendant le 1er confinement. 

addmargins(table(coco.activite$pcs, coco.activite$sitpro.conf, useNA="ifany"))
##                            
##                             En congé Au chômage En emploi <NA> Sum
##   Agriculteurs                     1          2         4    1   8
##   Artisans, comm., chef d'E        3          4        12    0  19
##   Cadres et pr.intell.sup.         5         14        59    1  79
##   Employés                        17         20        72    2 111
##   Prof. intermédiaires            21         26        55    3 105
##   Ouvriers                        10         20        38    1  69
##   <NA>                            21         51       124    3 199
##   Sum                             78        137       364   11 590
rprop(table(coco.activite$pcs, coco.activite$sitpro.conf, useNA="ifany"))
##                            
##                             En congé Au chômage En emploi <NA>  Total
##   Agriculteurs               12.5     25.0       50.0      12.5 100.0
##   Artisans, comm., chef d'E  15.8     21.1       63.2       0.0 100.0
##   Cadres et pr.intell.sup.    6.3     17.7       74.7       1.3 100.0
##   Employés                   15.3     18.0       64.9       1.8 100.0
##   Prof. intermédiaires       20.0     24.8       52.4       2.9 100.0
##   Ouvriers                   14.5     29.0       55.1       1.4 100.0
##   <NA>                       10.6     25.6       62.3       1.5 100.0
##   Ensemble                   13.2     23.2       61.7       1.9 100.0

Les effectifs sont faibles pour les agriculteur·rices et les artisans. On choisit de les exclure, ainsi que les enquêté·es n’ayant pas donné d’élément permettant de connaître leur situation professionnelle pendant le confinement.

coco.activite.khi <- coco.activite %>% 
  filter(! pcs %in% c("Agriculteurs",
                     "Artisans, comm., chef d'E")) %>% 
  filter(is.na(sitpro.conf)==F) %>% 
  filter(is.na(pcs)==F) %>% 
  mutate(pcs = fct_drop(pcs))
addmargins(table(coco.activite.khi$pcs, coco.activite.khi$sitpro.conf, useNA="ifany"))
##                           
##                            En congé Au chômage En emploi Sum
##   Cadres et pr.intell.sup.        5         14        59  78
##   Employés                       17         20        72 109
##   Prof. intermédiaires           21         26        55 102
##   Ouvriers                       10         20        38  68
##   Sum                            53         80       224 357
rprop(table(coco.activite.khi$pcs, coco.activite.khi$sitpro.conf, useNA="ifany"))
##                           
##                            En congé Au chômage En emploi Total
##   Cadres et pr.intell.sup.   6.4     17.9       75.6     100.0
##   Employés                  15.6     18.3       66.1     100.0
##   Prof. intermédiaires      20.6     25.5       53.9     100.0
##   Ouvriers                  14.7     29.4       55.9     100.0
##   Ensemble                  14.8     22.4       62.7     100.0

coco.activite.khi %>% 
  ggplot()+
  aes(fill=sitpro.conf, x = pcs)+
  geom_bar(position="fill")+
  labs(title = "Situation pro. pendant le confinement des personnes en activité avant selon la PCS",
       x="PCS",fill = "Situation professionnelle",y = "Proportion")

A première vue, les cadres et professions intellectuelles supérieures sont plus nombreux·ses à être resté·es en emploi, quand les professions intermédiaires ont davantage été mises en congés et les ouvriers et ouvrières plus souvent été au chômage pendant cette période. Il s’agit maintenant de vérifier ce résultat avec un test du khi-deux :

chisq.test(table(coco.activite.khi$pcs, coco.activite.khi$sitpro.conf))
## 
##  Pearson's Chi-squared test
## 
## data:  table(coco.activite.khi$pcs, coco.activite.khi$sitpro.conf)
## X-squared = 13.495, df = 6, p-value = 0.03582

Cette instruction donne le résultat pour l’ensemble du tableau :

  • X-square : la valeur calculée
  • df : le degré de liberté, soit la taille du tableau (nb de lignes - 1 x lb de colonnes - 1)
  • p-value : la probabilité que le lien observé soit du au hasard.

Par convention, les seuils suivants sont retenus :

  • p-value < 0,001 = très significatif, noté ***
  • p-value < 0,005 = assez significatif, noté **
  • p-value < 0,01 = un peu significatif, noté *
  • p-value >= 0,1 = non-significatif, noté °

Pour l’exemple, l’association observée est significative. On peut affirmer qu’il est probable que la pcs soit liée à la situation professionnelle durant le confinement des personnes qui étaient auparavant en activité.

…Une dernière chose :

Test du khi2 pondéré

# Produire le survey design :
coco.activite.khi.w <- svydesign(data=coco.activite.khi,
                                 ids=~1, weights = coco.activite.khi$POIDS)
# Tri croisé :
svytable(~ pcs + sitpro.conf, coco.activite.khi.w) %>% 
  rprop(2)
##                           sitpro.conf
## pcs                        En congé Au chômage En emploi Total 
##   Cadres et pr.intell.sup.  11.40    19.56      69.04    100.00
##   Employés                  14.29    17.29      68.42    100.00
##   Prof. intermédiaires      23.93    28.94      47.12    100.00
##   Ouvriers                  28.66    26.71      44.64    100.00
##   Ensemble                  20.23    23.74      56.03    100.00
svytable(~ pcs + sitpro.conf, coco.activite.khi.w) %>% 
  addmargins() %>% round(0)
##                           sitpro.conf
## pcs                        En congé Au chômage En emploi Sum
##   Cadres et pr.intell.sup.        7         12        43  62
##   Employés                       13         15        60  88
##   Prof. intermédiaires           28         33        54 115
##   Ouvriers                       22         21        34  77
##   Sum                            69         81       192 343
# Test du chi-2 pondéré :
svychisq(~ pcs + sitpro.conf, coco.activite.khi.w, statistic="Chisq")
## 
##  Pearson's X^2: Rao & Scott adjustment
## 
## data:  svychisq(~pcs + sitpro.conf, coco.activite.khi.w, statistic = "Chisq")
## X-squared = 19.516, df = 6, p-value = 0.2267

Cette fois-ci, le résultat tend à démontrer qu’il y a une forte probabilité pour que le lien visible entre les 2 variables soit du au hasard. Du fait des filtres que nous avons du faire, l’échantillon final sur lequel porte le test statistique est très faible.

Attention :

  • La logique du test du khi-2 est proche de celle des autres tests : on teste la probabilité de l’hypothèse 0 supposant une absence de lien entre les variables.
  • Les résultats sont dans les faits produits par les logiciels statistiques, il faut savoir les interpréter
  • Le test du khi-2 est symétrique -> même résultat quelque soit les modalités en lignes ou en colonnes
  • Le khi-deux dépend du nombre de modalités retenues -> on ne compare pas directement les valeurs du khi-2 mais les p-value
  • Si un effectif est inférieur à 5 -> opérer un regroupement de modalités et le test dépend du découpage de ces modalités
  • Le test est sensible aux effectifs -> à proportions égales, de plus grands effectifs seront plus significatifs
  • Le test (et la p-value) ne donne(nt) pas d’indication sur la force du lien entre les variables, uniquement si les variables sont probablement liées oui ou non -> voir le test du “V de cramer”
  • Le test du khi-2 ne donne pas de sens de causalité
  • En cas de p-value significative, on ne peut pas savoir si deux variables sont effectivement liées entre elles ou s’il y a une variable caché…

Plus de détails : Voir Barnier Julien, « Tout ce que vous n’avez jamais voulu savoir sur le χ² sans jamais avoir eu envie de le demander », 2016

Si le test du khi-2 est l’un des plus utilisés en sciences sociales (car il porte sur les fréquences), il existe une multitude de tests statistiques selon les configurations des données. Chacun porte ses propres conditions et hypothèses.

Voir aussi : https://larmarange.github.io/analyse-R/comparaisons-moyennes-et-proportions.html

Pour les estimations : https://larmarange.github.io/analyse-R/intervalles-de-confiance.html

Anova

Le test du Khi-Deux a une autre limite : il permet de tester l’indépendance entre deux variables nominales, mais comment fait-on pour mettre à l’épreuve la relation entre une variable qualitative et une variable quantitative (une variable nominale et une variable numérique) ? Parmi les tests possibles, le moins mal connu et le moins inutilisé en sciences sociales est l’ANOVA pour analyse de variance qui compare les variances entre différents groupes.

On compare la variance dans les groupes et entre ces groupes (intergroupe et intragroupe), les différents groupes correspondant aux modalités d’une variable catégorielle :

  • Si la variance intragroupe est supérieure à la variance intergroupe, on peut en conclure que le groupe d’appartenance ne modifie pas les valeurs de la variable numérique prise en compte : il n’y a pas de lien entre les variables ;

  • Si la variance intergroupe est supérieure à la variance intragroupe, alors c’est le signe d’une relation entre les deux variables ;

  • Sur le même modèle que le Khi-2, l’ANOVA compare les variances observées avec les variances théoriques en cas d’indépendance totale entre les deux variables ; l’hypothèse nulle est donc, comme pour le Khi-Deux, l’indépendance ;

  • Aux limites présentées pour le Khi-Deux s’en ajoute une autre : la variable numérique doit réunir certaines propriétés, notamment la normalité (avoir une distribution suivant une loi normale) ;

Dans R : anova <- aov(variableNumerique ~ variableNominale)

  • Où le ~ s’écrit AltGR+2 sous un PC, Alt enfoncé + N sous Mac) ;

  • Puis summary(anova)pour obtenir le résultat.

5. Production de rapport avec Rmarkdown

5.1. Qu’est-ce que Rmarkdown ?

Les fichiers et la synthaxe Rmarkdown permettent de générer des rapports automatiques mêlant du texte, du code et les résultats sur des jeux de données. Ce fichier contient :

  • Une en-tête définissant les propriétés générale du document produit (appelée YAML), soit les méta-données
  • Des titres permettant d’organiser le document,
  • Du texte brut, qui peut être mis en forme,
  • Du code (appelés chunks), exécuté ou non, présentant les retours ou non

Si vous êtes déjà dans un projet R (.Rproj), vous pouvez créer un fichier .Rmd :

Fichier > Nouveau fichier > Fichier R Markdown ; puis ok dans la fenêtre. Un fichier avec du texte d’exemple par défaut apparaît. Pour exécuter l’ensemble et produire le document : bouton “knit” ou ctrl + alt + K et nommer le fichier en sortie (uniquement la 1ère fois).

Exemple : Les présentations supports du cours ont été générées avec rmarkdown

5.2. Structure générale d’un Rmarkdown

En-tête (aussi appelée YAML)

L’en-tête est obligatoire et définit les principales caractéristiques du document produit. L’en-tête par défaut est minimale. Elle est donc définie en tout début de document et est délimitée par deux lignes “—”.

Attention, cet espace respecte l’indentation.

Elément nécessaire : Le format de sortie

output: html_document()

ou

output: pdf_document()

ou

output: word_document()

Voir : https://bookdown.org/yihui/rmarkdown/output-formats.html

Eléments fortement recommandés :

  • Le titre
title: "Titre du document"
  • Le(s) auteur(s)
author: "Nom de l'auteur"
  • Date
date: 'Dernière modification: `r format(Sys.Date(), format="%d/%m/%Y")`'

On peut voir que des instructions R apparaissent. Elles sont balisées comme ceci :

`r code R `

D’autres précisions utiles

  • Sous-titre :
subtitle: "Sous-titre"
  • Résumé :
abstract: "Ajout d'un résumé\nRetour à la ligne"

Des options propres au format de sortie :

Le format généré permet d’intégrer une série d’options dont :

number_sections: yes # pour ajouter des numéros aux titres
toc: yes # ajout d'une table des matières
toc_depth : '3' # Précision de la table des matières
folding_code : hide # Cacher le code par défaut
toc_float: true # Laisser la table des matières constamment visible

Et une série d’autres options précisées dans la cheat sheet Rmarkdown : https://raw.githubusercontent.com/rstudio/cheatsheets/master/rmarkdown-2.0.pdf

Exemple d’en-tête

---
title: "R - Avancé - M2 ENSP"
author: "Cécile Rodrigues - CNRS / CERAPS"
date: 'Dernière modification : `r format(Sys.Date(), format="%d/%m/%Y")`'
output:
  html_document:
    number_sections: yes
    toc: yes
    toc_depth: '3'
    toc_float: true
---

Rappel : L’indentation compte !

Ici, on a ajouté des numéros aux titres de sections, une table des matières allant jusqu’au 3ème niveau et présenté sur le côté (au lieu d’être en début de document).

Les titres

Les titres sont précédés par un # pour les titres de niveau 1 ; par ## pour les titres de niveau 2, etc.

Text brut

Le texte brut écrit dans le fichier .rmd apparaîtra comme le texte par défaut d’un traitement de texte. Le saut de paragraphe se fait en sautant deux lignes.

Les tirets seront matérialisés par des puces :

- Comme
- ici
    - En ajoutant des sous-catégories

1. Ou bien
2. avec des
3. numéros

La mise en forme italique est indiquée en entourant le texte à metttre en italique par * : *texte en italique* ; et celui en gras par ** (** texte en gras **)

Du code peut être inséré directement dans le texte :

`r 1+1 `

5.3. Les chunks (le code R ou autre)

Le code et les résultats R sont directement intégrés dans le texte mais balisé par

Les balises apparaissent en cliquant sur insert > R / python / etc ou en tapant Ctrl+Alt+I.

Exemple :

Sys.Date()

Donne :

## [1] "2022-11-23"

Plusieurs options peuvent être ajoutées dans l’en-tête du chunk :

{ r nom du chunk, eval = F, echo = F, warning=F }

Détail des options :

  • nom du chunk : Apparaît au même endroit que les titres de sections dans un script R classique
  • eval = F : n’exécute pas le code
  • echo = F : le code n’apparait pas
  • warning = F, message = F, error = F : les avertissement, message ou erreurs n’apparaissent pas
  • fig.width = 5, fig.height = 10 : taille des graphiques

Pour définir ces options systématiquement, ajouter en début de script .Rmd :

knitr::opts_chunk$set(warning=FALSE, eval=FALSE, message=FALSE, error=FALSE)

5.4. Tableaux et graphiques

Les tableaux et graphiques sont générés et mis en forme dans les chunks. Plusieurs packages permettent de produire des tableaux esthétiques, parfois complexes, et contenant des mises en forme particulières (couleur de fond, texte en gras, ajout d’une barre en fonction d’une valeur, fusionner des cellules, etc.). Dans le détail, ces fonctions produisent des tableaux en code html ou latex (suivant la fonction et le format de sortie) :

Les tableaux

un premier package utile pour faire des tableaux est {knitr} et les fonctions kable. Dans sa forme la plus basique, en reprenant l’exemple du tableau des “vies communes” :

vitavec <- freqm(coco, "q24", calc="%", exclude="6666", poids="POIDS") %>% 
  mutate(name= fct_recode(name,
                          "Conjoint·e"="q24_1",
                          "Enfants/petits-enfants"="q24_2",
                          "Parents/Grands-parents"="q24_3",
                          "Famille (autre)"="q24_4",
                          "Ami·es ou autres"="q24_5",
                          "Chien(s)"="q24_6",
                          "Autre(s) animal(aux)"="q24_7")) %>% 
  rename("Vit avec"="name", "Oui"="X1") %>% 
  select(-X2, -Total) %>% arrange(desc(Oui))
library(knitr)
vitavec %>% 
    kable("html", escape = F, caption = "Source : Coco1 - ELIPSS, CDSP, avril 2020")
Source : Coco1 - ELIPSS, CDSP, avril 2020
Vit avec Oui
Conjoint·e 81.6
Enfants/petits-enfants 48.5
Autre(s) animal(aux) 23.5
Chien(s) 14.2
Parents/Grands-parents 7.1
Famille (autre) 5.6
Ami·es ou autres 3.7

En un peu plus stylisé (en version html_document) :

library(kableExtra)
vitavec %>% 
    kable("html", escape = F, caption = "Source : Coco1 - ELIPSS, CDSP, avril 2020") %>% 
  kable_styling("striped", full_width = F)
Source : Coco1 - ELIPSS, CDSP, avril 2020
Vit avec Oui
Conjoint·e 81.6
Enfants/petits-enfants 48.5
Autre(s) animal(aux) 23.5
Chien(s) 14.2
Parents/Grands-parents 7.1
Famille (autre) 5.6
Ami·es ou autres 3.7

-> Améliorer les tableaux avec {kableExtra} : https://cran.r-project.org/web/packages/kableExtra/vignettes/awesome_table_in_html.html

Les graphiques de type ggplot deviennent des images mais vous pouvez agir sur leur taille. Ajouter les dimensions du graphique dans le chunk :

{r, fig.width= 9, fig.height=5}

vitavec %>% arrange(desc(Oui)) %>% mutate(`Vit avec` = fct_inorder(`Vit avec`)) %>% 
ggplot() + geom_histogram(aes(x = `Vit avec`, y = Oui), stat="identity", col="black",fill="white") +
  theme_bw()+ theme(axis.text.x = element_text(angle=45, hjust=1))+ 
  labs(y="Pourcentages (%)",
       title="Avec qui vivent les enquêté·es", caption="Source : Coco1 - ELIPSS, CDSP, avril 2020")

2 autres packages intéressants :

Les graphiques

Les graphiques de type ggplot deviennent des images mais vous pouvez agir sur leur taille. Ajouter les dimensions du graphique dans le chunk :

{r, fig.width= 9, fig.height=5}

vitavec %>% arrange(desc(Oui)) %>% mutate(`Vit avec` = fct_inorder(`Vit avec`)) %>% 
ggplot() + geom_histogram(aes(x = `Vit avec`, y = Oui), stat="identity", col="black",fill="white") +
  theme_bw()+ theme(axis.text.x = element_text(angle=45, hjust=1))+ 
  labs(y="Pourcentages (%)",
       title="Avec qui vivent les enquêté·es", caption="Source : Coco1 - ELIPSS, CDSP, avril 2020")

Les images

Enfin, les images peuvent être ajoutées de deux façons :

Dans le texte : ![Nom visible](chemin/vers/image/png)

Dans un chunk :

En ajoutant dans l’en-tête out.width = ‘400px’ pour préciser la taille :

knitr::include_graphics("img/tidydata_7.jpg")

5.5. Pour finir

Quarto

Cet été, la société qui produit RStudio a fait deux annonces. Elle s’appelle désormais posit mais le nom du logiciel reste rstudio et elle a lancé un nouveau projet : les fichiers quarto. Quarto est pensé comme un nouvel outil de publication scientifique et technique se rapprochant de rmarkown et des notebook (jupyter, par exemple). L’idée est de produire des documents hybrides, contenant du texte, de la rédaction, des commentaires et analyses mais aussi des figures, tables et résultats directement issus de codes intégrés dans le document. Quarto est présenté comme la “nouvelle génération” de rmarkdown et pousse l’outil sur d’autres langages (python, julia, javascript, etc.) ; ce qui était possible en rmarkdown sans que ça soit au centre de l’outil.

Voir : https://posit.co/blog/announcing-quarto-a-new-scientific-and-technical-publishing-system/

Un exemple : https://neocarto.github.io/docs/notebooks/bertinquarto/

–

Pour aller plus loin dans l’interactivité, vous pourrez utiliser le langage shiny qui permet de créer des applications web interactives dans l’éco-système de R. Cf l’intervention que vous aurez sur ce langage.

Raccourcis clavier

ctrl+alt+I : créer un chunk

ctrl+shift+K : “knit”, produire le document

Bonus. Exercices d’application

Comme on est passé·es un peu rapidement sur les exercices, prenons le temps de les faire maintenant :

Partie “Données pondérées” : Refaire les questions de rappel du M1 en version pondérée et correctement mise en forme :

  1. Commencez un fichier .rmd pour les exercices
  2. Quelle part de femmes et d’hommes ont répondu au questionnaire ? Présentez un tableau clair et lisible
  3. Faire le tableau et le graphique associé du sentiment de tristesse selon le genre des enquêté·es
  4. Quelle est la part moyenne des personnes de l’entourage des enquêté·es qui ont eu le covid et sont aujourd’hui guéries ? Et la part médiane ? -> En fonction de l’âge, En fonction de la PCS
  5. On veut maintenant comprendre l’effet du 1er confinement (et du début de la crise sanitaire) sur la situation d’emploi des enquêté·es
  6. Risques professionnels des personnes en activité selon le genre. Essayez de trouver une représentation graphique originale et adaptée à ce que vous voulez raconter !

6. Selon le temps, création de fonctions avec R

7. Création de fonctions avec R

Dans R, on crée des objets et on leur applique des fonctions pour produire des résultats et des analyses. Les fonctions prennent (1) un ou plusieurs arguments en entrée (ou paramètres ou input), éventuellement certaines options, (2) effectuent des opérations ou des actions et (3) produisent et renvoient des résultats ou sorties (output).

R est quand même un langage de programmation qui permet de créer ses propres fonctions si on ne les trouve pas dans le langage de base ou via un package.

L’écriture, dans sa version la plus basique, d’une fonction sera :

Ajoute2 <- function(x){
  res <- x + 2
  return(res)
}
Ajoute2(4)
  • x est l’argument, qui est donné en entrée par l’utilisateur·rice de la fonction
  • res <- x+2 est l’action : “Ajouter 2 au paramètre x”
  • return(res) renvoie le résultat, soit dans la console s’il n’y a rien à gauche, soit dans un objet si la fonction est précédée de objet <-.

Ajoute2 est apparu dans la partie “Functions” de l’environnement.

Attention : Si le nom de la fonction existait déjà, l’ancienne fonction est écrasée

summary <- function(x){
  res <- x+2
  return(res)}

summary est maintenant remplacée par le contenu de Ajoute2.

Mais on peut toujours supprimer la fonction summary qu’on vient de créer localement :

rm(summary)

Exemple d’une fonction qui renvoie le pourcentage d’un vecteur

propTable <- function(v){
  tri <- table(v)
  total <- length(v)
  tri <- round(tri/total*100,2)
  return(tri)
}
propTable(coco$q01)
## v
##     1     2 
## 19.89 80.11

Exemple d’une fonction qui renvoie un graphique

Ca peut être une fonction qui intègre la création d’un graphique :

MyBarplot <- function(var){
  tri <- table(var)
 barplot(tri, col ="skyblue", border=NA)
}
MyBarplot(coco$q01)

Pourquoi utiliser des fonctions ?

Les fonctions sont très utiles quand on se retrouve à répéter des actions. Par exemple, on recode une dizaine de variables sur le même modèle : 1 = Oui et 2 = Non :

Au lieu d’écrire…

library(questionr)
library(tidyverse)
freq(coco$q24_1)
coco$q24_1_r <- fct_collapse(factor(coco$q24_1),
                             "Oui"="1", "Non"="2", "NR"="6666")
freq(coco$q24_1_r)

… On crée la fonction qui fait ce recodage et on l’applique :

fct_ON <- function(x){
  recod <- fct_collapse(factor(x),
  "Oui"="1", "Non"="2", "NR"="6666")
  return(recod)
}
coco$q24_2_r <- fct_ON(coco$q24_2)
freq(coco$q24_2_r)

Les intérêts sont : - On évite une série de copier/coller qui alourdissent le code ; - Plus simple pour corriger d’éventuelles erreurs : on corrige la fonction et pas les 7 copies - Réutilisable d’un script à l’autre, vous pouvez faire des scripts avec vos fonctions habituelles

Comme pour les fonctions en général, les arguments peuvent être explicités ou être implicites selon leur position : fct_ON(x = coco$q24_2) ou fct_ON(coco$s24_2) car il n’y a pour l’instant qu’un seul argument.

En général, on compte sur l’implicite pour les données et on explicite les options; Comme dans mean(x, na.rm = T)

Un autre intérêt sera de tester un morceau de code en changeant un ou quelques paramètres :

graphede <- function(x){
des <- data.frame(jets=sample(1:6, x, rep=T))
ggplot(des)+aes(x=jets)+geom_bar(fill="white",col="black")+theme_bw()+scale_x_continuous(breaks = 1:6)
}
graphede(500)

Les valeurs par défaut

Il est possible d’intégrer des valeurs par défaut comme paramètres :

Nom <- function(x, dec=1, useNA="ifany"){}

–

Sans valeur par défaut, un argument est obligatoire mais facultatif si une valeur par défaut apparaît.

En ajoutant …, on indique qu’on intègre les arguments par défaut des fonctions appelées dans celle qu’on crée :

Nom <- function(x, ...){}

Objets locaux et objets globaux

  • Objet dans l’environnement = global / Objet dans la fonction = local
  • Une fonction peut accéder à un objet extérieur -> existant dans l’environnement
  • Mais les arguments et objets créés dans la fonction sont prioritaires
  • Un objet créé dans une fonction n’existe que dans la fonction
  • Les objets locaux sont détruits quand on sort de la fonction

Return

Par défaut, R renvoie le dernier objet présent dans la fonction. Mais en utilisant la fonction return, on peut davantage maîtriser ce qui est renvoyé. C’est par exemple utile quand on veut renvoyer un élément différent sur la base de conditions.

Comportement pas défaut :

fun <- function(x){
  x+2
  x+5
}
fun(2)

–

Pour forcer le retour de la 1ère ligne :

fun <- function(x){
  return(x+2)
  x+5}
fun(2)

Exercice

Créez une fonction qui renvoie le carré et le cube d’un nombre : ^ pour élever un chiffre à une puissance données : 2^2.

–

carrecube <- function(x) {
  carre = x^2
  cube = x^3
  return(list(x = x,
              carre = carre,
              cube=cube))
}

On renvoie une liste qui contient le nombre mis en entrée, son carré et son cube. Si on fait carrecube("Coucou") ?

–

On obtient une erreur qu’on peut éviter avec return car une fois qu’une fonction return est actionnée, le reste de la fonction n’est pas exécuté.

carrecube <- function(x) {
  if(!is.numeric(x)) return("Eh non. Il faut mettre un chiffre !")
     carre = x^2
     cube = x^3
     return(list(x = x,
                 carre = carre,
                 cube=cube))}
carrecube("coucou")
## [1] "Eh non. Il faut mettre un chiffre !"

Note : Si on met “print” au lieu de “return”, le code continue à s’exécuter et on rencontre une erreur.

Global / local

On n’est pas obligé·es de définir les arguments dans la fonction (local), ils peuvent être dans l’environnement global (la session rstudio). Mais si le même objet est dans la fonction, il sera priorisé.

Il est aussi possible de modifier un objet dans l’environnement global à partir d’une fonction en utilissant l’opérateur <<- :

x <- 4
changex <- function(a){x <<- a}
x
## [1] 4
changex(5)
x
## [1] 5

Exercice

Créez une fonction qui nous indique si on nombre est premier ou non

Un nombre est premier s’il n’est divisible que par lui-même et par 1, c’est-à-dire si le reste de la division est égal à 0 uniquement quand on le divise par lui-même ou 1.

  • Dans R, on connaît le reste d’une division à l’aide de l’opérateur %% : 4%%2
  • Une propriété des nombres premiers : Si un nombre n’est divisible que par 1 et aucun autre nombre premier inférieur ou égal à sa racine carré, il sera premier.
  1. Il faut envisager les arguments en entrée
  2. Il faut envisager les entrées possibles et donc les erreurs possibles -> Différents scénarios suivant ce que les utilisateur·rice mettent en entrée
  3. Ecrire les warnings en conséquence
  4. Envisager les cas de figure sans erreur et les sorties en conséquence

Créez une fonction qui renvoie tous les nombres premiers entre deux valeurs

premier <- function(x){
  # Pas tout seul-----
  if(length(x)>1) return(print("x doit être un nombre seul"))
  # Pas un nombre --------
  if(! class(x) %in% c("integer", "numeric")) return(print("x doit être un nombre"))
  # Pas entier ---------
  if(grepl("\\.",x)) return(print("x doit être un nombre entier"))
  # le nombre est-il premier ? -------
    for(i in (2:(x^.5))){
      if(x %% i == 0) return(paste0(x, " n'est pas premier"))}
    return(paste0(x, " est premier"))
}
premier(53)
## [1] "53 est premier"
listepremier <- function(x,y){
  liste <- NULL
for (i in (x:y)){
 for(j in (2:(i^.5))){
   if(i %% j == 0){liste <- c(liste, i)}
 }}
premiers <- setdiff(x:y, liste)
return(paste0(paste0("Les nombres premiers entre ",x," et ",y, " sont : "),
         paste(premiers, collapse = ", ")))  
}
listepremier(10,20)
## [1] "Les nombres premiers entre 10 et 20 sont : 11, 13, 17, 19"

Ressources

Point et aide pour le dossier

Rappel des consignes pour le dossier

Il s’agira d’un dossier en groupe de 2 à 4 étudiant·es au choix à partir de données fournies et à rechercher autour d’une même thématique. Le dossier rendu sera :

  • un document (.odt, .pdf ou .html)
  • rédigé de 4 à 6 pages
  • à partir d’un fichier .Rmd
  • comportant au moins une carte
  • présentant des résultats pondérés
  • et au moins un test statistique.

-> Vous m’enverrez le document produit et le script .Rmd au plus tard le 03/01/23 à mon adresse mail (vous aurez un accusé de réception sous 24h).

Attendus

  • Au moins 2 vagues de coco (conseil 1 et 6)
  • Résultats pondérés
  • Au moins un test statistique
  • Pas d’énumération de tests ou de tableaux sans réflexion, fil directeur et interprétation
  • Recherche de données géographiques pertinentes (voir data.gouv.fr)
  • Au moins une carte
  • Données transmises : coco - vague 1 à coco - vague 6, les autres vagues, fonds départements et communes

Discussion en groupes autour du sujet, d’un angle d’analyse à partir des données Coco (utiliser au moins 2 vagues).

Regarder la fiche “données covid” de data.gouv.fr

Rechercher les fonds de carte utiles

Ressources

R généralités

Tidyverse

Ggplot2

Cartographie

Aides en ligne et sites ressources

Quelques # et @ twitter “actifs”

Source : allison Horst @allison_horst

  • #rstats

  • #rstatsFR

  • @icymi_r

  • @hadleywickham

  • @thinkR_fr

  • @Rbloggers

  • @tidyversetweets